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ABSTRACT 



This thesis describes the study and development of a 
workable inverse scattering method for imaging and identifi- 
cation of radar targets. The space -time integral approach 
is used for iterative target shape reconstruction. 

Following an overview of transient electromagnetics, the 
integral equation is applied for thin-wire transient response 
computation. 

The analytical time domain integral equation is derived 
and solved numerically for general conducting bodies of 
revolution. Finally, the algorithm for an inverse scattering 
computer solution is derived and tested under simulation of 
physical environments. 
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I. 



INTRODUCTION 



A. OVERVIEW 

The classification of radar targets, has become increasing- 
ly important in the past two decades. In military applications, 
the need to distinguish between friendly, enemy, and neutral 
targets, is vital. 

In most radar applications, the only properties that are 
measured are range, angle, velocity and elevation. The methods 
that are used by radars to provide some classification of the 
target generally involve the examination of the echo signal. 

Some of these techniques can be summarized briefly as follows: 

1. High range resolution radar, which uses a very short 
pulse that can provide sufficient range resolution to obtain 
a rough profile of the illuminated part of the target shape. 

2. Cross section fluctuations analysis which can provide 
information on the size of the target. 

3. Synthetic aperture radar that gives a high angular 
resolution. This method is typically employed in ground 
mapping . 

4. Analysis of the backscattered polarization from a target 
can provide a means for discrimination between types of 
targets . 

5. Use of nonlinear effects due to junction of metals. 

The backscattered spectrum is analyzed and checked for 
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harmonics appearing in the received signal, that are typical 
for certain type of targets. 

6. Doppler modulation by propellers, the rotary wing in 
helicopters or even hull vibrations in armored vehicles. 

None of the methods cited can be used to obtain all the 
information on the target shape and size, they only can be 
used to discriminate to a certain extent between different 
types of targets. Other methods currently used involve 
passive reception of radiation by targets. 

Elint receivers or IR detectors can be used to monitor 
radiation from targets, and, based on previous intelligence, 
or signature catalog, associate a type of radiation to a type 
of target. 

The inverse scattering technique , that is the subject of 
study in this thesis, uses the basic properties of the transi- 
ent response to an impulse excitation (impulse response) in 
order to obtain information on the size, shape and also material 
composition of the target. 

In principle, the smoothed impulse response can be ob- 
tained by measuring the backscattered fields at all frequen- 
cies up to some upper limit. In most cases this would be 
impractical so, instead, the target is illuminated by a very 
short pulse (with a very wide band frequency spectrum ranging 
from the low frequencies in the Rayleigh region up to higher 
resonance region frequencies) . 
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The impulse response contains all the information needed 
to characterize the target. The higher frequencies give the 
information about the structure. This is related to the 
optical theory of diffraction which is a good approximation, 
under certain conditions, for shape prediction up to the 
shadow boundary. Physcial optics, however, fails to predict 
the phenomena appearing in the shadow region, such as the 
creeping wave which is due to lower frequency waves travelling 
around the body. 

A different approach that has been used to retrieve tar- 
get configuration information is related to the singularity 
expansion method (S.E.M.). This method seeks to classify 
scatterers through the complex residues and poles of their time 
domain signatures. Briefly, in this method the transient E.M. 
response of a conducting body is characterized as a series of 
complex exponentials (or damped sinusoids) : [1] , [2] 

N St 

f(t) = I e 

m=l 

where R is the amplitude (residue) of each mode of complex 
m 

frequency s^^. The true function has a complex frequency 
representation of the form: 

N R 

, . r ITl 

“ ^^s-s 

m=l m 

The poles s depend only on the geometry (independently 
of the orientation) , and the residue depends on both 
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excitation and geometry [3]. Thus^once the poles and residues 
are known, they provide a means for target classification. 

The problem is to find the poles and residues. 

A method is suggested using Prony's algorithm in [2] for 
extraction of poles and residues from time-domain signatures. 
The method has been tried successfully on wire geometries and 
simple shapes. 

Another important application of the impulse response of 
a target is the study of its radar cross section (RCS) , since 
the Fourier transform of the impulse response, (the frequency 
response) is directly related to the RCS by the following 
relationship : 

a = IT c^ ]f (w) I ^ 

where 

F(0)) = / F^(t) e^“^ dt 

»oo 

and 

Fj(t) is the impulse response of the target. 

Finally the time-domain analysis of the impulse response of 
a body is of great importance in the study of the electro- 
magnetic pulse (E.M.P.) response of targets. 
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B. PRACTICAL CONSIDERATIONS 



The practical study and measurement of the impulse re- 
sponse involves the building of a time-domain scattering range. 
The availability of generators having very short pulse-widths, 
fast rise time, and high voltage outputs, together with fast 
sampling oscilloscopes and small computers makes the task 
feasible . 

For the solution of the inverse scattering problem the 
time domain scattering range developed at the Naval Post- 
graduate School [4] , was employed to measure the impulse 
response and demonstrate the validity of the theories developed 
in this thesis. The transient response data measured by using 
pulse generators with very fast rise times can validate the 
calculation of the time domain integral equation solution. 

More important, it provides a vehicle for testing various tar- 
get imagery and I.D. schemes, as well as investigating broad- 
band RCS reduction. 

C. SURVEY OF THE LITERATURE 

One of the first efforts to determine physical properties 
of electromagnetic scatterers from time-domain analysis was 
due to Kennaugh and Cosgriff [5] , who employed the physical 
optics approximation to calculate the approximate backscattered 
impulse response of a flat spheroid. 

A summary of this work was reported by Kennaugh and Moffatt 
[6] , who extended the physical optics approximation to transient 
response calculations. At this stage the computation of the 
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impulse response was approximated by using the physical optics 
approach or the inverse transform of the frequency domain solu- 
tion. Then, in 1968, Bennett derived in his Ph.D. thesis [7] 
an approach for solving the time domain integral equation 
directly by using a time -stepping method. Bennett used a 
form of the equation called the magnetic field integral equation 
(or M.F.I.E.) for two and three-dimensional surfaces. 

Sayre and Harrington [8] used the same time-stepping tech- 
nique to derive the transient response of a wire, but they 
employed another type of equation called the electric field 
integral equation (EFIE) . 

Additional work in time-domain electromagnetics has been 
reported in several more recent papers as, for example, the 
report by Miller, Poggio and Biirke [9] on time-domain analysis 
of thin-wires. Summaries and review of the time-domain electro- 
magnetics are due to Mittra [2] and [10] and Miller [11]. 

A comprehensive review on the subject and its application 
has been published by Bennett [12] , which gives the reader 
an introduction and initiation to time-domain electromagnetics. 
Literature dealing with the inverse scattering problem is due 
to Young [13] and Shubert, Young and Moffatt [14], which des- 
cribes the image generation of a target from the ramp response 
waveform. 

Moffatt and Mains [15] describe a method of detection and 
discrimination of radar targets using the idea of ramp response 
but with multiple harmonic radar frequencies. Finally, a very 
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recent report by Bennett [16] gives the general approach to 
the time domain inverse scattering solution. The method sug- 
gested in that report is being used in this thesis to approach 
the problem of transient target imaging. This list of efforts 
is by no means exhaustive, as many other authors have published 
papers dealing with transient electromagnetics and target 
imaging . 

D. THESIS OBJECTIVES 

»The objective of this effort is to develop the analytical 
and numerical solution of time-domain inverse scattering 
problem. It follows the development of the time-domain 
scattering range at the NFS [4]. 

The main effort has been conducted toward the formulation 
and numerical solution of the time-domain integral equations 
for conducting bodies of revolution. 

A computer program is developed to solve for direct and 
inverse scattering. The algorithm was used for experimental 
display of the shape of the target under test. 
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II. GENERAL SOLUTION OF THE TIME 
DOMAIN INTEGRAL EQUATIONS 



This chapter is a review of the development of time domain 
integral equations , their solutions and the associated numeri- 
cal considerations. Detailed development of the equations 
can be found in [2,10,11]. 

A. THE INTEGRAL EQUATIONS (E.F.I.E. and M.F.I.E.) 

The basis for derivation of the equations are the general 
time dependent forms of Maxwell's equations : 



V • E(r,t) = p(r,t)/e^ 

V • H^(r,t) = 0 

Starting from these equations, and using the expressions for 
the magnetic vector potential. A, and the electric scalar 
potential , c|) , which are related to the fields by 



V xE(r,t) 




V xH(r,t) 



It- E(r,t) + J(r,t) 

at O 






- 7 (()(?, t) - ||(r,t) 
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one can obtain the expression of the magnetic field integral 
equation (M.F.I.E.) for a perfect conductor [2] 



J3(r,t) 



/N ->T_nr' ^ 

2n X (r,t) 



^ ^4 dS’ <2-l) 



R 



and the electric field integral equation (EFIR) 



n X E 






R 9x s 



2— <1-^ 



q (r ' , t) (1 , 1 9 • 

£ R^ 
o 

In equations (2.1) and (2.2), 



(2.2) 



n is the unit vector normal to the conductor at 
the observation point r. 

J (r',x) is the current density at the source point r' 
at a retarded time x . 



q(r',x) is the charge density at source point at the 
retarded time x . 



Here x is given by 



X = t - R/c 



and 



R is the distance from the observation point r to 
the source point f ' , i.e., 



R 



r 



r 



I 
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E (r,t) and H (r,t) are the incident electric and magnetic 
fields at the observation point. Because of numerical con- 
straints, the equation most suitable for closed surfaces is 
the MFIE equation (2.1) while for thin wires and open struc- 
tures the EFIE equation (2.2) is used. 

Equation (2.1) states that the current induced on a point 

“^inc 

r of a body at a time t is due to the incident field H at 
that point and the sum of scattered fields there from earlier, 
more distant locations of the body. The interactions between 
the samples of the body are displaced in time by an amount 
equal to the transit time of the fields between them. 

In the integral, the observation point (R = 0) is excluded 
since this "self-term" is included in the incident field por- 
tion of the equation. This determines the general solution 
for the MFIE equation which is solved directly by marching 
on in time. The solution is started at a certain time (t = 0) , 
and due to causality the integral part of the equation will 
be zero since there are no past currents. Then, as t increases 
the integral is solved with previous currents already computed. 
This technique is a time-stepping procedure for solving initial 
value problems . 

For the EFIE, equation (2.2), the method is slightly differ 
ent but uses the same time stepping procedure. The solution 
of this equation is applied to the problem of transient 
scattering by a thin wire, in the next chapter. 

The solution of equations (2.1) and (2.2) gives the cur- 
rents induced on the body by the incident field. Once the 
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currents are known, it is a straightforward step to compute 
the scattered fields. For the far-field case the field is 
given by: 



;>scatt ->• 
H 





(2.3) 



where a_ is the unit vector from the source r ' to the far 



observation point and r^ is the separation distance. J^(r*,x) 
is the current at the source point at a previous (retarded) 
time T = t - r /c. 



B. NUMERICAL SOLUTION AND CONSIDERATIONS 

The solution of the integral equations implies finding 
the currents induced on the body once illuminated by an inci- 
dent field. Since the solution is to be carried out numerically, 
the first step is to represent the currents by their space and 
time samples J^.(s,t) in a discrete manner. J^j(s,t) is the 
current at the i segment AS (or patch) at the j time 
increment AT. The total current will be approximated in terms 
of all samples at all times and spatial segments via the basis 
function expansion: 



where Ng is the number of space samples on the body while N^ 
is the total number of time increments, and the pulse basis 
functions are: 



o 



J(r,t) = I I J..(s,t) U(t.) V(s.) (2.4) 

i=l j=l ^ 
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1 at the center of the jth time interval 



U(tj) 



0 elsewhere 



1 at the center of the ith space patch 



V(si) 



0 elsewhere. 



The space time sampling is shown in figure 2.1. The 
sampling of the currents already suggests that the body is 
divided into segments (or patches) AS and the time is sampled 
at AT intervals, while the total observation time goes up to 

X at. 

It is very important at this point, before attempting to 
solve this equation, to define the basic criteria for proper 
selection of the space and time sampling as well as their 
relationship to the shape of the exciting pulse and the body 
size . 

The incident field is chosen here to be represented by a 
gaussian impulse due to its rapidly decreasing features and 
also because most of the short pulse generators will produce 
such a waveform. 

The gaussian impulse is given by 



g(t) 




e 



Its amplitude falls to 1/10 of its maximum value at 
tmax “ 1.5/A. Its frequency spectrum is given by [2] 



G(f) = F{g(t)} 



e 
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Figure 2.1. Space-tirae sampling diagram 
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This response reduces to 1/10 of its maximum value at 
^max “ 0*5A. The time sampling AT must exceed the Shannon 
rate to avoid aliasing problems if FFT's are employed later, 
i .e . , 



AT < 
min 



1 



2f 



max 



1 

A 



Usually AT will be taken smaller. 



AT 



5f 



max 



1 

2.5A 



(2.5) 



The spatial sampling AS must be chosen such that 



AS ^ ATxc (2.6) 

otherwise there will be coupling between the currents of 
adjacent patches. Violation of (2.6) will also mean that the 
stepping procedure for solving the equations, which assumes 
that during one time interval the currents at one patch are 
not influenced by adjacent currents, will no longer be satis- 
fied, thus producing an incorrect solution. 

Another consideration is that the space sample points must 
be close enough to resolve adequately the spatial variation 
of the incident field over the body. From Eq. (2.5) and (2.6) 
we have: 



AS 



- 5f 



max 
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AS > 



A . 
min 

5 

This latter relationship imposes some minimum ratio of the 

size L of the body to the minimum wavelength (or equivalently 

the minimum pulse width) . The ratio L/A . determines the 

min 

resolution and the usefulness of the time domain analysis, 
and it must be kept high enough to maintain adequate resolution. 
If the body is to be represented with good fidelity with a 
certain minimum number of samples N we will need 

Ne A . 

L > 

— 5 

This means that if A . increases (i.e., wider pulse) L must 
increase or equivalently for a small size target A^^^ must 
decrease (i.e., shorter pulse). 

Having presented the equations, their general solution and 
the requirements for the various parameters involved we can 
move on to the implementation. 

C. APPLICATION TO THE TIME -DOMAIN SCATTERING BY A THIN-WIRE 
The first application that was studied was the computation 
of the impulse response of a straight thin wire. This example 
was chosen since it gives a good understanding of the physical 
processes and the numerical techniques used in the solution 
of time domain integral equations. The basis for this compu- 
tation is given by Sayre and Harrington in [8]. 

The problem is formulated as follows: a thin-wire of 

length L and radius a aligned with the z axis is illuminated 
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by an impulse field of the form 



ginC(^) _ ^ max 



The basic formulas are in terms of the retarded potentials; 
The boundary condition : 



E^'"‘"(t,z) 



3A 



z , 3<j) 



8t 3z 



where the retarded magnetic potential is given by 



A^ (t,z) 



= JL f ,T) ^ , 

4n ■> P 



Axis 



R 



and the retarded electric potential is 



<P = 



4TT£ 



o Axis 



f d 

'* R 



where the continuity equation between charge and current is 
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3z 



- p. 

3t 



In these equations I(z',x) and q(z',x) are the current and 
charge at source point z' at time x = t - R/c , while R is 
the distance between the source and observation points. The 
geometry is described in figure 2.2. 

Using the numerical formulation expressed in [8], and 
applying the time stepping procedure, a computer program was 
written (Appendix D) to compute the scattered field from a 
thin wire of length 1 meter and radius 0.8 mm. The wire was 
divided into 40 segments and the sampling time AT = AZ/c = L/40c. 



24 




Figure 2.2. Geometry of the thin-wire scatterer 
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The incident field was taken to be similar to the pulse 
transmitted by the generator of the NPS transient scattering 
range, with linear polarization parallel to the wire and 
incidence normally on it. 



ginc 




(t-t 



max 



2 



A 



6 . 10 ^ 



-1 

sec 



t = 18 'AT (about half transit time along 

““ the wire) 



The program computes the currents at each spatial segment with 
up to 400 time samples. This corresponds to 10 transit times 
across the wire. The current at the center of the wire is 
given in figure 2.3. The backscattered field is given in 
figure 2.4. 

No significant change was obtained for 20 segments and 
200 time increments. Comparison of this computation with lab 
measurement was done by Hammond [4] and very good agreement 
was found . 

Similarly, comparison of the computed results found in 
[11] with a ratio of L/a = 0.00667 and A = 3.25*10^ gave 
exact agreement. The backscattered field for this case is 
shown in figure 2.5. 
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CURRENT (MILLIRMP5) 




TIME (LIGHT-METER] 



Figure 2.3. Current induced at the center of a 

thin-wire by a gaussian incident pulse 
(A = 6- 108 and a/L = 0.0008) 
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NORMRLIZED FIELD (MILLIVOLTS) 



400 . 




-400. ‘■ 

0 . 



_J I I \ I I I I L_ 

2. 4. 6. 8. 10. 



TIME (LIGHT-METER) 



12 . 



Figure 2.4. 



Backscattered field from a thin-wire 
(A = 6-10^, a/L = 0.0008) 
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NORMRLIZED FIELD (MILLIVOLTS) 



lOO 



<2 [11] Figure 58 




TIME (LIGHT-METER] 



Figure 2.5. Backscattered field from a thin-wire 
(A = 3.25.1 q 9 and a/L = 0.0067) 
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III. TIME DOMAIN RESPONSE OF PERFECTLY 
CONDUCTING BODY OF REVOLUTION 



A. INTRODUCTION 

This section considers the solution of the time domain 
response for a body of revolution, when the exciting wave is 
propagated along the axis of symmetry of the body. 

The use of the axial symmetry of the body will result in 
considerable simplification of the computation of both the 
induced currents on the body surface and the field scattered 
from the body. The algorithm developed in this section will 
be used later for the solution of the inverse scattering 
problem. 

B. THE MAGNETIC FIELD INTEGRAL EQUATION 

In [7] Bennett derives the integro-dif ferential equation 
for three dimensional scatterers: 




2a xH (r,t) 
n 




where : 



T 



t - R/C (retarded time) 



r is the observation point 



r' is the source (or integration) point 
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R = I r - r ' 1 

A is the surface of the scatterer 

H^(r,t) is the incident magnetic field at the observer’s 
point. 

The geometry for the body of revolution is given in figure 

3 . 1 . 

For the actual problem the incident magnetic field is 
propagating along the z-axis (axis of rotation of the body ) , 
and is directed towards the positive x-axis. The body lies 
completely to the right of the positive z-axis, and is defined 
by the continuous contour function, p(z), which is the radius 
of the body at a point z. 

The prime coordinates denote the source or integration 
point, while the unprimed denote the observation point for 
which the current is desired. At each point on the body, there 

/\ /\ /N 

is defined the triad of orthogonal unit vectors, 
where : 

/\ 

a : unit vector normal to the surface 

n 

/N 

a . ; unit vector tangent to the surface and directed 
towards the positive 4) direction 

a^; unit vector tangent to the surface and directed 
along the curve representing the longitudinal 
contour of the body. 

For the following derivation we will define the angle a (or 

/N /N 

a') as that sustained by a (or a') with the z axis, a (or a') 

s s 

being positive when a (or a') points outward from the z axis 

S 5 
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and negative if it points towards it. The angle (j) (or (j) ' ) 
is the angle defined by p (or p') with the (x,y) plane. 
Equation (3.1) states that the current density at the obser- 
vation point r at time t is due to: 

1. The current induced by the incident field at (r,t) . 

This part represents the physical optics approximation, 

2a xH^(r,t). 
n ' ' ' 

2. The current due to the flow of currents at all other 
points of the body at a retardated time ,t = t - R/C, where 
R is the distance from the integration point to the oberva- 
tion point. 

According to the geometry described in figure 3.1 we can write 



The current on the body will flow only in the s and <)) direc- 
tions and can be expressed (for both primed and unprimed 
coordinates) , as 



the current components and are functions of the position on 
the scatterer. We will denote them simply by J and J, so 



R 



(p^ + p'^ - 2pp' cos (())' - ())) +(z' - z)^) 



J (r ,t) 



Jg(p,z,t,(j))*ag t J^(P/Z/t,(j))»a^ 




s 



that 



J(r,t) 
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Figure 3.1. Geometry of the scattering by 
a body of revolution 
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Eq, (3.1) can be written in terms of its components in the 
s and (j) direction, i.e., 

J (r,t) = [2a xH^(r,t)]. i |-] [2 x j' x r]^ dA* 

S n g Z7T K d T n S 

(3.2) 

= [2a^xgi(r,t)Jj + + xg).dA' 

(3.3) 

In order to obtain the scalar expressions of the two last 
equations we have to perform the vector multiplications; 

This is done in Appendix A and as a result, equations (3.2) 
and (3.3) become ; 



= -2HNz,t)sin4. -j'v-singJdA' 

A R ^ 

(3.4) 

= -2H^(z,t)cos4.sina+i /4[i|p + il[J'Usin 8 +j;F 2 ldA- 

A R 

(3.5) 



where ; 



S = (J,' - (j, 

= (z' -z) cos 3 sin d + (p -p' cos 3) cos a' 

V = z' - z 

= (p cos 3 “ p' ) cos a + sin a cos 3 (z ' “ z) 

U = (z ' - z) sin a • sin a' +p cos a sin a ' - p ' sin a *cos a ' 
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We note that the incident field is dependent only on z and 
t (independnet of 4* due to synunetry) . We also note that; 

dA' = p ’ d<ti ' ds ’ = p’ d3 ds’ 

The integrals in Eq. (3.4) and (3.5) can then be divided into 
two parts; integration over B from 0 to 2tt, and integration 
over ds' from 0 to S (S being the contour curve as shown in 
figure 3.2). 




The next step is to obtain further simplication of the inte- 
gral equations (3.4) and (3.5). 

The currents in these equations are functions of their 
position on the body; 

dg — jg(Z,P,(f),t) 

= J^(z,p,4-,t) 

But we remark that the currents due to the incident field have 
a sinusoidal variation with >p , as follows; 

= H^(z,t) cos ♦ sin a 
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= H^(z,t) sin <i> 

o 



so that we can assume that the currents will have the same 
respective relationships with (|). 

J (z,p,t,(j>) = J (z,p,t)sin (J) = J *sin (j) (3.6) 

s s s 

J, (z,p , t,(j)) = J.(z,p,t)cos (J) = J,*cos <p (3.7) 

9 q> 9 

This assixmption is demonstrated in Appendix B, using Maxwell's 

equation. 

After substitution of equations (3.6) and (3.7) into the 
integrals (3.4) and (3.5), the integrals over the angle B can 
be simplified as shown in Appendix C, so that the scalar 
integral equations reduce to their final expressions; 






[J'F, COS 3 + Vsin^BJl ] d3 
Si (p 



(3.8) 



S 2tt 

= -2H^(z,t)sinc+ i p'ds' / 

0 R 



[J' U sin^3 + cos 3 J;]d3 

S 2 (p 



(3.9) 



in which the currents J and J, (J' and J') are only functions 

s 9 s 9 

of p, z, and t (p',z',t) and independent of the angle (j) . 

The significance here is that the surface distribution of 
the unknown currents is reduced to a one-dimensional variation 
of current along the generating arc of the surface of revolution. 
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This will allow considerable reduction in computer time and 
storage when carrying out the numerical solution of the 
equations . 

Once the currents along the contour of the body are known, 
the currents at any other point (at some time t and some z 
position) will be given by 

Jg(PfZ,t,(J)) = Jg(p,z,t) -sin (p 
(p ,z,t,(f)) = (p,z,t) -cos 4>. 



C. FAR -ZONE SCATTERED FIELD 

In most cases we are interested in the far scattered 
field. At a large distance from the body the scattered field 
is given by: 



H®(r,t) 



4irr C 
o 



A 






X Sj^ldA' 



(3.8) 



where 



T is the retardated time t - R'/C. 

R' is the path difference between the integration 
point to the observer and the origin to the 
observer. 

J(r',x) is the current at the surface integration point 
retarded in time, while a^^ is the unit vector 
pointing from the source point towards the obser- 
vation point. The geometry is described in 
figure 3.3. 



r is the distance from the origin to the observation 
° point. 
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A is the projection of the source point on the 

plane (x,z) 



Figure 3.3. Geometry for the far-zone 
scattered field 
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In the following, the scattered field will be derived 
for an arbitrary elevation angle 9. Due to symmetry with 
(J) there is no loss of generality if we observe the scattered 
field in the plane (x,z) where (}) = 0 . 

According to the geometry: 

R* = ~[z' cos 9 + p' cos (J) ' sin 9] 

T = t - Pf/C. 

R' is negative whenever the observation point is closer to the 
source than the origin, and positive otherwise. 

/N 

a„ = sin 9 a + cos 9 a 

K. X Z 

Performing the vector multiplication in equation (3 . 8 ) and 
separating the field into its cartesian components, gives the 
followoing set of three scalar equations ; 

H®(9,t) = -! — 7 :r / [J sin^())'sin a' +J , cos^(}) ' ] cos 9 dA' 

X 4 Trir^^ oT s (p 

O A ^ 

(3.9) 

H®(0,t) = / |:p[ (J, - Jg sina')sin (})' cos (})' cos 9 

^ o A ^ 

+ J sin (j) ' sin 9 cosa'JdA' (3.10) 

s 

H^(9,t) = - 5 — j I — [- sin 9 ( J_sin^(J) ' sin a' +J , cos^cj) ' ] dA ' 

Z 47T]r C oT S CD 

O A ^ 

(3.11) 

After substitution of dA' = p ' d (}>' ds ' , and using the odd 
properties of some of the terms appearing in the integral on 
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(p ' (as it has been done previously in Appendix C) , the three 
scalar equations representing the components of the scattered 
field reduce to the following expressions; 



nj(9,t) 


= cos 9 H®(t) 


(3.12) 


X 


0 


(9 ,t) 


= - sin 9 H®(t) 

0 


(3.13) 


Hy(9,t) 


= 0 


(3.14) 



3 

where the common term H^(t) is given by the following integral 

“o'*"’ ' J J J 3 ' sin a’ 

+ cos^4)')d4)' (3.13) 

Equations (3.12) to (3.15) represent the scattered field for 
an arbitrary elevation angle 9 . The special case of interest 
is the backscattered field where 9 = tt . For this case there 
is only one component of the scattered field, which is 

Hj(t) = - H®(t) 

X o 

and the retarded time is: 

T = t - R'/C 



where 



R' = 
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All the terras appearing inside the integral in equation 
(3.15) are independent of (}) ' so that the integration on cj) ' 
can be carried out analytically. Doing this yields the 
following expression for the backscattered field: 



H (t) = 

s 



4r C 

o 0 



/ (p ' ,z ' ,t) sin a' +J^ (p ' ,z ' ,t) ] p 'ds ' 



(3.16) 



In equation (3.16) , the currents J and J, are those given 

S (p 

by equations (3.8) and (3.9) respectively, used at tirae 
X = t - z'/C (z* ^ 0) , where z'/C is the retardation tirae rela- 
tive to the origin. The backscattered field equation is 
independent of (p which raeans that we can integrate the tirae 
derivative of the currents, along a line on the contour of the 
body. 



D. NEAR ZONE BACKSCATTERED FIELD 

In the Tirae Doraain Scattering Range developed at NFS, 
the distance from the antenna to the target under test is not 
long enough to be considered for all frequencies of interest 
as a far region. Consequently, it is needed to develop a 
forraulation of the scattered field in the near zone. For this 
case the georaetry is described in figure 3.4. 

The general equation for the scattered field is given by: 

H®(J,t) = ^ xS]a|i (3.17) 
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point 



Figure 3.4. Geometry of the near-zone 
backscattered field 
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where 



X 


= t - R'/ C , and 


dA' 


in 

Q. 

II 




R 2 = (r^ + z'^) + p'^ 

O / K 



After carrying out the vectorial multiplication, dividing the 
field vector into its cartesian components , again dropping 
all the odd terms in (j) ' and integrating analytically on (f) ' , 
equation (3.17) reduces to; 



H®(r ,t) = 

X o' 


4 J „,2‘r'-"c 3^HJ3(P' cosa’ - (r^ + z')sin «• 

U R 

- J , (r + z ' ) ] P 'ds ' 

* ° (3.18) 



Equation (3.18) is the expression of the near backscattered 



field where; 






R' = { (r + z ') ^ + P 

o 



The currents J {p',z',T) and J, (p',z',T) are given by equations 

S ^ 

(3.8) and (3.9), being used at the retard time t. 

X = t- (R'-r^)/C (referred to the origin) 
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IV. NUMERICAL SOLUTION OF THE TII4E DOMAIN 
RESPONSE OF THE BODY OF REVOLUTION 



In the previous section the analytical equations for the 
currents and scattered field were derived. The equation ob- 
tained will have to be solved numerically, and the first step 
will be to divide the body into patches of approximately equal 
curvilinear surface area. This is followed by the setup of 
the algorithm for the solution. 



A. DIVIDING THE BODY INTO PATCHES 

The geometry of the body is completely specified by the 
function p(z), (0 < z ^ L) ; which represents the body contour 
function S. Since the integration is to be carried over S, 
the contour will be divided into NS segments of equal length 
AS. Consequently, the body contour function is fully defined 



by NS+1 


points . 








n 


The space samples of 
segment and is given 


each : 
by 




"i = 


^*^£+1 




z . = 

1 


(^£+1 



£ = 1 to NS+1 



£ = i = 1 , . . . ,NS 



(4.1) 



+ z ^)/2 



The angle a for each sample point is the slope of the curve 
at that point and is approximated linearly by 
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Pp+l 

a . = Arc tan { ) , £ = i = l,...,NS (4.2) 

^ ^l+l ~ 



Figure 4.1 shows the geometry of the contour sampling. Each 
segment AS defines a ring at with radius p^. All the rings 
will be divided into patches of approximately equal surface 
area, dA^, according to the following formulas* 

dA. = DS • p. • Ag. = AS^ 

1 11 

where A3^^ is the angle between two adjacent sample points on 
the same ring. 




This suggests a division of each ring into an even number of 
patches NPI which will be 

TTp . 

NPI = 2xINT(-^) (4.3) 



The function INT means the integer value of its argument so 

*t h 

that the increment angle for the i^ ring becomes 




2it 

NPI 



(4.4) 



Thus, following this procedure, any arbitrarily shaped body 
of revolution can be subdivided into patches of approximately 
equal surface area. Figure 4.2 shows the geometry of the 
space sampling on the whole body. 
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Figure 4.1. Geometry of the contour sampling 
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Figure 4.2. Space sampling of the whole body 
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B. NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS 

The equations to be solved are given by Eq. (3.8) and 
(3.9) in previous sections. 

For the numerical solution the currents will be represented 
by their space-time samples: 

J(P/Z^t) — iT(i/n) 

where i is the i space sample on the contour function and 
n is the n^^ time sample 

i = 1 , . . . , NS 

n = 1 , . . . , NT 

t = n-AT 

The time being sampled at a rate AT. Similarly the incident 
field is sampled and will be denoted by its n-th time sample 
at the i-th segment: 

H^(z,t) = H^(i,n) 

The general procedure for the numerical procedure was already 

mentioned in Chapter II, and will be reviewed in more detail 

here. The incident field, H, is assiomed to be zero at all 

times less than some reference value (which will be set equal 

tolAt). Hence the total current on the surface of the scatterer 

is also zero prior to time t^. 

At the next time interval t, = t + At the incident field 

1 o 

reaches the scatterer at z = Zj^ (first space sample) , thus at 
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time the current will be simply twice the tangential com- 
ponent of the incident field. As time marches on, the incident 
field covers more of the body and the current will be due to 
the incident field plus the fields due to currents at other 
points on the scatterer at earlier times. So, by marching 
on in time, the currents at each space sample, for each time 
increment, can be computed. 

The currents appearing inside the integrals in equations 
(3.8) and (3.9) are the retarded currents , and are given by 

The point (p',z') is represented by the k-th space sample 
while the time t will be: 

T = t - R/C = n-At - R/C = gAT. 

In the general case x will not be an integer multiple of 
AT, so that the current cannot be represented by any of the 
previously computed samples. In order to find the value of 
the current at time t , there is a need to interpolate between 
the sampled known, discrete values. 

For a better accuracy use was made of the four point 
Lagrangian interpolation formula using two samples at each 
side of the value gAT if t , t^r t^ , t^ are the four samples 
used for interpolation, such that: 

t^ = t^ + At = + 2At = + 3At 

^4 = g4 At 
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At 



= gj At 
= g^ At 

g^/ g^r are integer values, and 
T = g At (g real number) 



such that 



ti < t2 < T < t3 < t^ 



The interpolation formula using those four points will be 
given by: [17] 



4 T-t. 



S , (p 



I I ” (r— 1 

n=l i=l Vi " 

iT^n 



(4.5) 



This interpolation formula must be adjusted in two cases: 

1. Whenever t^ or t^ is greater than the actual time 
of computation t (in order not to interpolate in 
the future) . 

2. Whenever t^ or t 2 are less than the beginning time 
of computation t^ (due to causality) . 

The time derivative of the currents will be simply expressed 
as the derivative with respect to x of the interpolated 
value, which will give the following formula: 
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(4.6) 



s,(p 

3 t 



4 T-t, 



I Irl n (i — . (t ) 1 

nil i=l Vi " 

i^n 



Defining; 



Tl = g - 9 i 
Tj = 9-92 
Tj = 9-93 
T4 = 9-94 

The interpolation formulas for the current and its time 
derivative can be written in their numerical form as follows 






and 



DJ 



S,(J) 



T X T X T- T_ X T- X T 

P J .(ifg.) - ;= 

6 s,(j) ^4 6 



T, X T . X 

2 '^s,*<^' 92 ) 



T4 X Tl X T^ 






(4.7) 






3 x 



T„T_ -l-T-,T. +T„T, 



6 AT 






6 AT 




T 3 T 1 +T 4 T 1 +T 3 T 4 


T 


2AT 


^S,(j) 


T.T, + T, T_ + T_T . 




4 1 12 2 4 


.T 


2AT 


^s,(j) 



(4.8) 
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Now we can write the numerical expression for the currents 
in equations (3.8) and (3.9) by representing the continuous 
integration by niomerical summation, and expressing the cur- 
rents in terms of their space-time samples . 

For simplicity, we will write first the expression for the 
integral part of equation (3.8) and (3.9) which will be 
denoted by J , and J . . 

CS C<P 



NP. 



J^s(k,n) 



, NS i 

= ^ J, I 

1=1 m=l 



mji^l for I=k kim 



V*sl^ (m) J (i,g) +Fj^CO(m) Jg (i,g) 



kIm 



V-sI^ (m)DJ. (i,g)+F,CO(m)DJ_(i,g) 

+ i ^ s 1 



(4.9) 



NP. 






, Mb 1 

i I Ae.PjASj I 
1=1 m=l 



mT^l for i=k 



R^. 

kim 



O.SI (m) (i,g)+F2C0(m) J (i,g) 

• I - - — ■ - - ^ 

R, 



kim 



U- SI^ (m) DJg (i,g) +F2CO (m) DJ^ (i ,g) 
+ _ 



] 



(4.10) 



where ; 



J (i,g) are given by equation (4.7) and DJ (i,g) 
by equation (4.8). 



kim 



(p^+p^-2pj^p^CO(m) + (z^ -Zj^)^)^/^ 



CO(m) = Cos [ (m-1) A3j^] 
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SI (m) 



= sin [ (m-1) A6^] 

V = Z. - z 
X u 

U = V sin sin + Pj^ cos sin ot -p^ sina ^ cosa 
= V*CO (m) • sin a . + (p, -p . CO (m) ) cos a. 

jl 3. JC X X 

= COS ( Pj^cos (m) ” Pj_) + V sin a^^CO (m) 

NS represents the total number of segments 

NP . represent the total number of patches at the 
ith ring. 



The total currents can now be expressed in terms of J and 



and are the values of the current components stored and 
used for computation of the backscattered far field, which 
will be expressed in numerical form as follows; 



Jc^f which yields; 





(4.11) 




-2H^(Zj,,n) X sin (a^^) +J^^(k,n) 



(4.12) 




NS 



H^(n) 



o i=l 



where 



g 




and DJ^ ^(i,g) are given by equation (4.8). 
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Equations (4.7) and (4.8) used for interpolation of the 
currents and their time derivative are given for the general 
case where; 

1. g^^ > 1 (after the incident field had reached the 
specific point i on the scatterer) 

2. g^ < n (so that no interpolation in the future occurs) . 

In cases where these two conditions are not met the interpo- 
lation must be adjusted to include only those points which 
satisfy causality with no interpolation in the future. 

The numerical formulations derived in this section were 
programmed in FORTRAN to compute the currents and the back- 
scattered field, due to the excitation of a body of revolution 
by an impulse type field (or a ramp) . The input data to the 
program is the function p(z) along with the length of the 
spatial segment AS and the time increment AT. The program 
is given in Appendix E. 

C . NUMERICAL EXAMPLES 

The first example to be considered is a sphere of radius 
1 meter. This special example was chosen since its result 
appears in various publications and comparison can hence be 
made . 

For the present example , use is made of the guidelines set 
in Section (LI .B) for the choice of the incident field, the 
sampling rate and the space sampling . The incident field for 
the case of the impulse response is represented by the gaussian 
shape 
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H-"(z,t) = 



( t-t -z/c ) ^ 

^ max 2a 

e , — < t < 0 



, otherwise 



where 



A = 6*10^ sec ^ 



which yields 



AT = 0.2/C sec 



The pulse is described in figure 4.3. 




Figure 4.3. Gaussian shaped pulse 
incident on a sphere 

The pulse reaches its maximum value at the center of the sphere 



t 

max 



a/C . 
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The spatial segment, with 15 segments 



AS = = 0.209 m < C*AT 

The number of time increments observed was such that it 
covers 6 times the transit time across the whole sphere (60-AT). 
The results obtained for the backscattered field are given in 
figure 4.4, on which the results of Bennett [7] are superimposed 
for comparison. 

Observation of the impulse response of the sphere shows 
the following results: 

1. The first peak appears after a transit time across 
one radius of the sphere (which is the time where the inci- 
dent pulse reaches its peak value) , and is due to the return 
from the leading edge of the sphere . 

2. The second peak appears at about 5.2 a/C and is due to 
the return from the back side of the sphere, after a transit 
time of : 



Si 

c 



+ 




a 

C 



5.2 



a 

C 



which corresponds to twice the free space transit time up to 
the shadow boundary plus the transit time around the back 
of the sphere. This second peak is the creeping wave return. 

3. The response as expected, starts decaying after some 
time equivalent to 5-6 times the transit time across the sphere. 

The ramp response for the same sphere is given in figure 
4.5. The first positive incursion of the response waveform 
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BRCKSCRTTERED FIELD (RRNGEmH) 
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Figure 4.4. 



Smoothed impulse response of a 
1 meter radius sphere 
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BRCKSCRTTERED FIELD (RRNGExHxC) 
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Figure 4.5. Ramp response of a 1 meter 
radius sphere 
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extends to about 5 times the transit time across the radius 
which corresponds to 1.25 round trips across the target. The 
impulse response of a 30 cm diameter sphere was measured in 
the NPS scattering range lab. The waveform obtained is shown 
in figure 4.6. The response appears here with a reversed 
sign since vertical polarization is used (instead of the 
horizontal polarization considered for the computations) . 

The waveform shows very good agreement with the form obtained 
for the 2 meter diameter sphere. 

The second example considered was the cone-sphere whose 
geometry is shown in figure 4.7. The target contour was 
divided into 12 segments while all other parameters remained 
the same as for the sphere . The impulse response of the cone- 
sphere target is given in figure 4.8. The waveforms show two 
positive peaks, corresponding to the leading edge and creeping 
wave returns. The separation between them is equivalent to the 
two way free space transit time up to the shadow boundary 
plus the transit time across the rear of the sphere. 
[(2xl.3+7rx0.7) = 4.7 light-meters] The ramp response 

of the same cone-sphere target is given in figure 4.9, and 
the remarks made on the sphere apply here too . 

The third example was the capped cylinder target shown in 

figure 4.10. The contour of the target was divided into 12 

segments, the sampling time was AT = 0.025 m/sec and 
9 -1 

A = 4*10 sec . The smoothed impulse (gaussian) response 
of this target is given in figure 4.11. The waveform shows 



59 



clearly the responses from the leading edge and the rear of 
the cylinder, with the same amplitudes, while from the flat 
horizontal region of the cylinder the response is zero as 
expected from the analytical expressions of the field. 

The ramp response for the capped cylinder is shown in 
figure 4.12. The first positive incursion of the response 
clearly exhibits a shape comparable to the shape of the target. 
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Figure 4.6. Measured impulse response of a 30 cm diameter 
sphere (incident field vertically polarized) 



Radius 




Figure 4.7. Geometry of the cone-sphere 
target (half-target) 
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Figure 4.8. Smoothed impulse response of 
the cone-sphere 
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Figure 4.9. Ramp response of the cone-sphere 
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Figure 4.10. Geometry of the capped cylinder 
target (half -target) 



bHLKbCHTTERED FTELD (RflNGExH) 



0. OM 



0.02 



0,00 



- 0.02 



-0.04 




L -L-J-I. I I I I I. . I I I I I. I I I I I I I I I 1 I I ,,, I , 

0.00 0.25 0.50 0.75 1.00 1.25 1.50 




TIME (LIGHT-METER) 



Figure 4.11. Smoothed impulse response of a 

capped cylinder (length ~ 0.33 m) 
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Figure 4.12. Ramp response of the capped 

cylinder (total length ~ 0.33 m) 
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V. INVERSE SCATTERING 



A. PROBLEM DEFINITION 

Given the incident field and the scattered field from an 
unknown target, it is desired in the most general case to 
retrieve the shape, size, orientation, and, perhaps, composi- 
tion of the target. This is known as inverse scattering, or 
target imaging. 

For the present case, the problem is restricted to a 
metallic body of revolution illuminated by a field propagating 
along the axis of symmetry, and the scattered field is measured 
in the backscattered direction, as shown in figure 5.1. 

The solution of this problem is approached here directly 
in the time-domain. 



B. FORMULATION OF THE PROBLEM 

In [6] , Kennaugh and Moffatt derive the relationship between 
the physical-optics impulse response of a body and its area 
function A(z) , projected in the plane orthogonal to the direc- 
tion of propagation of the field (z-axis) . 

The impulse response is given by: 



Fj (t) 



8irr 



d^A(z) 

dz^ 



z=ct/2 



(5.1) 



where r^ is the range from the observation point to the 
scatterer, as shown in figure 5.1. Integrating twice (5.1) 
with respect to time gives the ramp response as follows: 
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Figure 5.1. Inverse scattering geometry 



69 





z=ct/2 



(5.2) 



Equations (5.1) and (5.2) were proven to be valid for bodies 
having a monotonic area function with a limiting value at the 
shadow boundary. This relationship between the ramp response 
and A(z) is derived in Appendix F for the body of revolution. 

The total backscattered field that would be measured will 
consist of the physical optics term plus an additional term 
due to the field induced interactions between the currents 
at different parts of the body. The total current induced on 
the body is given by equation (3.1). The first term of the 
R.H.S. of the equation represents the physical-optics current, 
which we will denote by J^, and the second term expressed by 
the integral represents the additional current, which we will 
denote by J^. The total current will then be: 



Expressed in terms of (5.3) the total backscattered field for 
a ramp excitation will be: 



J 




5 (r,t) + J„(r,t) 



P 



(5.3) 




1 




(5.4) 



47rr C 
o 



H®(t) + H®(t) 
P c 



where : 



H^(t) is the physical optics ramp response 
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3 

H^(t) is the additional response, 

and 

H„(t) is the total ramp response. 

Using the relationship between the physical-optics response 
and the area function, (5.4) can be rewritten as: 

H^(t) = 2 ^ , z = ct/2 (5.5) 

o 

which leads to 

p^(z) = 2r^C(H|(t) - H®(t)) (5.6) 

Equation (5.6) is the basis for the solution of the inverse 
scattering problem for the body of revolution . 

C. INVERSE SCATTERING SOLUTION PROCEDURE 

Equation (5.6) relates the contour function of the body, 
p(z), to its ramp response. Solution of that equation leads 
to the shape of the body. The only term known in equation 
(5.6) is the measured value of the total ramp response Hj^(t) 
while (t) is not known. 

Moreover, all the terms in the equation are functions of 
p(z), hence there results a non-linear equation that will be 
solved by an iterative procedure. As a first approximation, 
the unknown term H^(t) will be neglected and so the first 
estimate of the contour function will be given by: 

p^(z) = [2r^C H|(t)]^/^ (5.7) 
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Once a body is defined its ramp response can be computed since 

we know the incident field. The ramp response computed for 

the first estimate of the body will be denoted by H (t) and 

^1 

is expressed as: 

2/ X 

s s 

< 5 . 8 ) 

1 O 1 

A new estimate of the shape 92^"^^ found using (5.6) 

and (5.7), and is given by 



2/ X 

P2(z) 



2r C[H;(t) - 
O R 



(t)] + 



2/ X 

(z) 



(5.9) 



or equivelently by 



p^iz) = 2r^C[2H®(t) - H|^(t)] (5.10) 

The new estimate P 2 (^) serves to compute a new ramp response, 
(z) , from which an updated estimate of p(z) is extracted. 

^2 

Following this iterative procedure, the k-th estimate of the 
contour function will be given by: 



P^(z) = Pk_i(z) + 2r - h| (t)] (5.11) 

k-1 



or by: 



2 / X 

p„(z) 



k-1 

2r C[k.H®(t) - I Hp (t)] 
° ^ £=0 



(5.12) 



The iterations are continued up to the point where the differ- 
ence between the true measured response and the last computed 
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response reaches a minimum. Several runs of the program for 
various bodies has shown that, as the iterations go on, the 
mean squared error is decreasing, reaches a minimum then starts 
oscillating. In one case observed, the error is diverging 
after reaching the minimum. Some representative results are 
shown in the next section. 

The normalized mean squared error will be defined as 
follows ; 



e 



2 

k 



(H®{t) - H® (t))^ 

__2 \-l 

(H®{t) 



(5.13) 



2 

e, is the value that must be minimized in order to obtain the 
k 

most accurate representation of the shape of the body. As 
it is defined, the error is equivalent to the difference be- 
tween the successive shapes, as shown by equation (5.11). Once 

the difference H (t) - H„ (t) is minimized, the difference 

^ "^k-1 

2 2 ^ 

Pj^(z) - P]^_]^(z) is also minimized. 

D. NUMERICAL SOLUTION AND EXAMPLES 



For the numerical solution we will use equation (5.12), 

which is more suitable for the algorithm developed. 

We will incorporate the following terminology; 

Hj^j^(nAT) is the time sample of the measured ramp response 

while (nAT) is the time sample of the computed ramp response 
^k 

for the k-th estimate of the body, and P]^(z) is the k-th 
estimate of the contour function. In discrete numerical 
form, equation (5.12) becomes: 
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(5.14) 







k-1 

I H (nAT)] 
£=0 



where 



z 



n 



n C AT/2 



and for £ = 0, H (nAT) = 0. The mean squared error is: 

rv^ 



N 



2 




I[H„^(nAT) - H (nAT)] 
n=l ^ ^k-1 



(5.15) 




where N is the total number of time samples observed and 
computed. 

The ramp response of the body defined by will 

be computed using the integral equation solution to the for- 
ward scattering problem, as already described in Chapter IV. 

The Fortran program, given in Appendix E, has been extended 
with a subroutine named "SHAPE", in order to implement the 
inverse scattering algorithm. In this subroutine the shape 
of the body is computed iteratively, using (5.14). 

The space samples given by (5.14) are defined at equal 
increments on the z-axis, while the program developed for the 
computation of the ramp response has been proven to work best 
when the contour was divided into equal length segments. Con- 
sequently, one of the tasks of subroutine "shape" is to find 
the space-samples at the points delimiting segments of 
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approximately equal length. This is done using a "cut and 
try" method which converges to the desired length of segment. 

The iterative approach for solving the inversion equation 
(5.14) was tested in a simulative way for the three types of 
bodies already considered for the forward scattering solution 
in the previous chapter. 

The ramp response of the body is first computed and is 
used to simulate the response that would have been measured 
(without noise or error) . Then the shape of the body is re- 
trieved from that ramp response. 

The shape of the body is retrieved from the first positive 
incursion of the ramp response, which as it will be seen 
extends beyond the length of the target for the shapes 
considered . 

The first example was the 1 meter radius sphere . Figure 
5.2 shows the 1-st, 3-rd, 6-th estimates of the contour func- 
tion of the sphere. The 6-th estimate gives the least error 
and is shown superimposed on the original shape in figure 5.3. 

The second example was the cone-sphere, used in two con- 
figurations. The first one when the contour of the body was 
originally cut to 12 segments. Figure 5.4 shows the successive 
estimates (1,3,6) of the contour, and figure 5.5 shows the 
minimum error estimate (estimate 6) . The second version, 
when the body was cut to 11 segments , gave a much more accurate 
result, as shown in figures 5.6 and 5.7. 
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Figure 5.2. Contour estimates of the sphere 
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Figure 5.3. 6-th estimate of the contour of the sphere 
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Figure 5.4. Contour estimates of a cone-sphere 
(cut to 12 segments) 
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Figure 5.6. Estimates of the contour of the cone 
sphere (oust to 11 segments) 
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Figure 5.7. Sixth estimate of the contour of the 
cone sphere (cut to 11 segments) 



The last example was the capped-cylinder . Figure 5.8 
shows the successive estimates of the shape (estimates 1,3,4), 
and figure 5.9 shows the 4-th estimate which gives the minimum 
error. 

For the case of the capped cylinder, after reaching the 
minimum error at the 4-th iteration, the imaged shaped of 
the body started to diverge considerably, and appeared more 
shortened than the actual shape . 

Common to all examples, the error appears mostly at the 
shadowed region of the body, but in general the shape is 
retrieved with quite a good fidelity, which proves that the 
iterative direct time domain solution is working, subject to 
certain caveats which will be discussed in the conclusions. 

E. INFLUENCE OF NOISE 

In the examples considered in the previous section the 
only noise appearing in the ramp response is the numerical 
noise generated by sampling and computational round-off. In 
this section the measurement noise will be tested with a ramp 
response to which noise is added. In this simulation, each 
noise sample is obtained by the sum of a large number of zero- 
mean, uniformly distributed random numbers, which, by use of 
the central limit theorem, approximate the statistics of 
gaussian noise. The noise samples are normalized in order to 
remain within a certain voltage signal to noise ratio, and 
then are added to the computed ramp response. 

The body considered for the noise test was the sphere, 
and the signal to noise ratios tested were: 
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-th approximation 
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V 

20 log — = 10, 20, 30 and 40 dB . 

^noise 

The sphere ramp responses with the noise added are shown in 
figures 5.10, 5.11, 5.12 and 5.13 for the successive signal 
to noise ratios. 

For the 10 dB S/N case, the appearance of spikes and dis- 
continuities on the ramp response causes the shape of the 
body to diverge considerably since the first iteration, and 
the shape cannot be recovered. 

For the 20 dB S/N case, the computed shape starts to con- 
verge , with the minimum error shape obtained at the third 
iteration. Then the imaged shape proceeds to diverge and to 
shorten . 

Figure 5.14 shows the successive estimates of the shape 
up to the 3-rd, and figure 5.15 shows the sixth estimate, which 
is clearly different of the actual shape. 

For the 30 dB S/N case, the minimum error shape was ob- 
tained at the 5-th iteration, as shown by figure 5.16. Then 
again, the computed shape diverges as shown in figure 5.17. 
Common to all cases of divergence, the wrong estimates start 
as soon as the last computed shape exhibits a strong discon- 
tinuity in its contour. A very important conclusion from the 
noise analysis is that before starting the inverse scattering 
procedure the ramp response used must be smoothed so that no 
spikes or sudden discontinuities appear on it. 

We should normally expect that the true measured response 
will have a 20 dB S/N ratio, so that the problem of divergence 
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Figure 5.10. Ramp response of the sphere with 
noise added (S/n = 10 dB) 
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Figure 5.11. Ramp response of the sphere with 
noise added (S' /N = 20 dB) 
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Figure 5.12. Ramp response of the sphere with 
noise added (S/N = 30 dB) 
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Figure 5.13. Ramp response of the sphere with 
noise added (S/N = 40 dB) 
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Figure 5.15. Sixth estimate of the shape of the 
sphere (S/N = 20 dB) 
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Figure 5.16. Estimate of the shape of the sphere 
(S/N = 30 dB) 
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Figure 5.17. Eighth estimate of the shape of the 
sphere (S/N = 30 dB) 



must be investigated. In any case, before diverging the 
shape that gives the minimum error is quite a good approxi- 
mation of the true shape. 

An additional test with 40 dB S/N was made, and no signi- 
ficant deviation, from the case without noise added, was 
observed. Figure 5.18 shows the first and 6-th estimates of 
the shape. The minimxam error shape was the 6-th estimate. 

Up to the 12-th iteration the shape of the body remained 
stable, with minor variations. 
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Figure 5,18. Estimates of the shape of the sphere 
(S/N = 40 dB) 



VI . SUMMARY 



A. SUMMARY OF RESULTS 

The inverse scattering method for radar target imaging 
was the main topic addressed in this work. The inverse 
scattering algorithm developed uses the time domain solution 
of the integral equation. The class of targets considered 
included the conducting bodies of revolution whose important 
symmetry property reduces considerably the computation time 
and the computer storage needed for the solution of the in- 
verse scattering problem. 

The use of the ramp response for retrieving the shape of 
the target proved to be a very efficient tool, and the exam- 
ples considered gave generally good results, with initial 
convergence of the imaged shape. For the case of the long 
cylinder and the sphere with low signal to noise ratios there 
appeared an eventual divergence of the imaged shape after 
many iterations. This is thought to be due to accumulated 
numerical roundoff errors in repeatedly solving the integral 
equations. In general it has been shown that the use of 
the ramp response by itself, without any further computation, 
gives a close approximation to the shape of the body but up 
to its shadow boundary which corresponds to the maximum of the 
area function A(z) along the line of sight. Beyond that boundary 
the body is, at the first observation, elongated and the 
iterative algorithm is needed to compensate for the errors 
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induced by the wake or "creeping wave" which alters the 
targe t re spon se . 

As the iterative computation process goes on, some numeri- 
cal instability that is mainly due to a discontinuity in the 
computed shape might arise. However, before diverging, the 
shape of the body always reaches an estimate which has a 
minimum error relative to the true shape, and is a good 
approximation of that shape, so that the minimum error cri- 
teria set for stopping the iterative process appear to be an 
attractive way to overcome the instability phenomena. As a 
continuation of this work it is suggested to carry out a 
detailed analysis of the error and numerical instabilities 
that have been observed while testing the iterative algorithm. 

B . RECOMMENDATIONS 

The algorithm developed for target imaging by transient 
inverse scattering, has been tested numerically by simulating 
additive physical noise. It must now be evaluated in a real 
laboratory environment where the true measured response is 
used. The transient scattering range at NFS should be used 
to test and finalize the algorithm for the class of targets 
addressed in this thesis. 

The first step to be done is to obtain the measured ramp 
response of the target. The range uses an impulse generator 
which exhibits a gaussian pulse shape characteristic. The 
true ramp response can be obtained by using either a time 
domain convolution or a frequency domain computation and I.F.T. 
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The true ramp response can be obtained, in principle, from 
the gaussian response of the target using the following 
re la tion ship : 



Hr(03) 



Hig(03) 

G(w) 



2 



where : 



Hr(w) is the Fourier transform of the true ramp 
response of the target, 

H__(w) is the gaussian impulse response of the 
target, 

G(w) is the Fourier transform of the transmitted 
gaussian impulse, 

and the ramp response will be: 



HR(t) = f"^[Hr(w)] . 

Very special care must be given to the choice of space 
sampling and time sampling intervals. The time sampling is 
dictated by the sampling rate of the oscilloscope, and the 
space sampling must be done at a high enough rate such that in 
every case the distance R between the closest spatial samples 
is related to the time sampling T by: 

R/C ^ AT . 

The next step is to translate the FORTRAN program given in 
Appendix E to fit the computer of the lab. There is a need 



98 



to reduce to a minimum the computation time since the 
laboratory microcomputer is very slow compared to a large 
mainframe. The use of some additional symmetries of the 
body can certainly help in this task. 

Finally, it is necessary to minimize the noise in measure- 
ments, and smooth the response to be used for shape computation. 

C. AREAS FOR FURTHER STUDY 

The ultimate objective of this research program is to 
aid in determining the feasibility of developing future radar 
systems having the innate ability to classify targets using 
their transient scattering response. 

In the most general case the target will be of arbitrary 
shape and the use of the iterative method described in this 
thesis would be highly impractical and prohibitive for com- 
plex targets, due to two principal factors. The first is the 
time needed to compute the response of a complex body, and 
the second is that the area function obtained from the ramp 
response is not sufficient to determine in what way the 
response is to be computed. 

The best technique to deal with such targets would involve 
a way that does not require the computation of the integral 
equation, and instead will use only the measured impulse 
response (or ramp response) to discriminate between targets. 

The first solution involves measuring the ramp response at 
multiple look angles. The response waveform will be used up 
to the shadow boundary to find the area function along each 
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line of sight. The target image can be approximately recon- 
structed from the projections of the area functions. Das 
and Boerner [18], suggested a very interesting approach for 
target shape estimation using these projections. This method 
seems attractive and is felt that it could lead to a good 
solution for target imaging. Its main disadvantage is the 
need for multiple look angles, requiring either interfaced 
multiple radars, or long-time dwelling with one radar. 

An alternative approach is to classify targets via the 
complex poles and residues of their time-domain signatures. 

Once the poles and residues are extracted they could be com- 
pared with stored library data corresponding to known targets, 
and thus would provide a classification of the target londer 
test. This catalogue approach is also being pursued using 
the NFS transient range, as well as through other independent 
efforts, and may prove to be the most practically implementable 
technique available. 

Another area which could form the subject of future re- 
search includes the time correlation of the shape of the 
target with its motion. 
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APPENDIX A 



DERIVATION OF THE SCALAR INTEGRAL EQUATIONS 

In Chapter III, equations (3 . 3) and (3 . 4) involve some vector 
multiplications. The first appearing in the integrals is: 



->• ->• 
a ^ j ^ R 
n 



(a xa' xR).j' + (a x a ' xr)j' 
n s s n (j) (j) 



(Al) 



We note the following identities ; 



a X a ' X R 
n s 



xs xs 



a' (a .R) 
s n 



- R(a .a') 
n s 



(A2) 



a X a ' X R 
n (j) 



a ' (a .R) 

cp n 



R(a •a') 
n (j) 



(A3) 



^ ^ ^ /N /N 

Expressing R, a' and a' in terms of the (a ,a.,a ) coordinates 

s 9 s cp n 

system gives; 



^ /N /X 

R = a, (R*a.) + a (R*a ) 

9 (p S S 

S XS ^ ^ 

^ x\ ^ ^ ^ ^ 

i' = a,(a'*a ) + a (a'*a ) 

s 9 s s s s s 



(A4) 



For the body of revolutions we note the following relationships 



XN XX 



9 9 

^ XV 


cos 3 


a • sin 


B 


* a = 
9 s 


- sin 


^ . XV 


a! *a = 

9 n 


- COS 


a • sin 


3 




a' -a 

s s 


sin a 


sin a 


'cos 3 + cos a cos a' 
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sin a' sin 3 



a ' • a 

s <J) 

/s ^ 

a ' *a 

s n 



= cos a sin a 'cos 3 - sin a cos a' 



where 



3 = (p • - (p 



The vectors R, a , a , a, are given in cartesian coordinates 

n s (|) ^ 



by: 



R = a ( p cos <}) - p ' cos <J) ' ) + a ( p sin <t> - p ’ sin (p ' ) 
X y 



- a (z ' - z) 
z 



i = a (cos a *cos d)) + a (cos a *sin d)) - a sin a 

n X y ' ^ z 

/N /N 

I = a (sin a *cos <J>) + a (sin a *sin <J> ) + a cos a 

s X y ' ^ ' z 



= -a (sin a) + a *cos d) 

X y 



Carrying out the dot products gives 



R*a 



n 



= cos a(p -p'cos 3) + (z' -z)sin a 



R*a = sin a(p - p'cos 3) - (z' -z)cos a 
s 



(A5) 



R* a . = - p ' sin 3 



Using the relationships given by the sets of equations (A4) 
and (A5) , and simplifying, equations (A2) and (A3) become: 
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/V /N /N 

a xa' xR = a. {sin g [ (z ' -z) sinasina ' +pcosasina' 

ns <p K I V / K 

- p 'cosa ' sina] } + a { (p-p 'cosg) cosa ' 

s 

+ (z ' -z) cosgsina ' } (A6) 



/\ /N ->• 

a X a ! X R 
n (j) 



/N 

a^ [ (z '-z) sinacos g+ cosa (pcosg-p ' ) ] 

+ a [-sing (z ' -z) ] 
s 



(A7) 



Defining : 



V = z ' - z 

= (z ' - z) cos g sin a ' + cos a ' (p-p 'cos g ) 

U = (z ' - z) sin a sin a ' + p cos a sin a ' - p ' cos a ' sin a 

= cos a(pcos g -p ' ) + sin a cos g(z' - z) 



produces from (A6) and (A7) 

a xa' xR = a. U sin g + a *F, (A8) 

ns (|) s 1 

a X a ' X R = a. F- - a V*sin g (A9) 

n q) <p 2 s 

Note that V, Fj^, F 2 are even functions of g and U is inde- 
pendent of g . Next we have to carry the second vector multi- 
plication : 



a X = -a sin (|) - a. cos (p sin a (AlO) 

ns <p 
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where is given in cartesian coordinates by: 

After substituting equations (A8) , (A9) , and (AlO) into 

equations (3.3) and (3.4) we obtain the scalar expressions 
for the components of the current in the s and ()) directions 
as follows : 



J = -2HSin4, + i + vsln eldA’ (All) 

A R ^ 

_ Ay.r I > , 1 r 1 r 1 . 1 3 . 

J* = - 2 H cos ♦sin a +27 /t-'r-'c 

A R 

• [J1 U sin 3 + J!F„]dA' (A12) 

S (p 2 
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APPENDIX B 



DEMONSTRATION OF THE SINUSOIDAL VARIATION 
OF THE CURRENTS WITH j) 



Maxwell's first equation; 



V xH = e E + J 



(Bl) 



In the case of interest, due to syiranetry, the field H is 
independent of (j) and is given by 



H = a H (p,z,t) 

X o 



or in cylindrical coordinates 



H = H cos d) a - H sin cb a, 
o P o 4> 



8H 



V X H = [— -^{p H ) ]sin (j) a^ + *cos ()> a, 

pSz'^o ^ p 3z ^ (p 

• , o 

- a sin (b ’-s; — 
z 3p 



(B2) 



The currents on the body expressed in cylindrical coordinates 
are : 



y\ 

J = a sin a J + a, J, + a cos a J (B3) 

P S (|> (|) z s 

The electric field; 



E = a E = a sin (b E + a, cos cb E (B4) 

yo p ^ 0 (() ^o 

Si±>stituting equations (B2) , ;(B3) and (B4) in Maxwell's 
equation (Bl) and equating the terms in the same direction gives 
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1 9 

sin a'J = [— — (pH ) - E ] sin <|> 

s p o z o o 

The term inside the brackets is independent of <}> so that we 
can write 



sin a J = J (p,z,t) sin d) 
s °i 

Similarly 

9H 

- E )cos <}> 

<}> 9z o ^ 



J = J (p,z,t) cos (p 

Y ® O 



and finally 

9H 

J cos a = - sin a 

s 9p^ 



(B5) 



(B6) 



= J (p,z,t) sin <}> (B7) 

°3 

Equations (B5) , (B6) , (B7) show clearly that the currents 

Jg and vary sinusoidally with the angle <}> so that we can 
generalize : 



J (p,z,t,(}>) = J (p,z,t) sin 4> 

S o 

(p ,z,t,(}>) = J^(p,z,t) cos (p 
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APPENDIX C 



SIMPLIFICATION OF THE SCALAR INTEGRAL EQUATIONS 

In Chapter III the scalar equations for the currents 
components were given by equations (3.4) and (3.5). In these 
equations we had: 

dA' = p' ds' d6 



and 



Js = J sin (f) 
s 

Jcj) = J. cos (f) 

<f> 

Substituting in equation (3.4) gives 

sin ♦ = -2H sin ♦ + 5 ^ J p'ds' J + g-I 

• [J ' sin (})'F, - J ' sin 6 V cos 4> ' 3 (Cl) 

S J. (p 

As noted in Appendix A, the expressions V, F^ and F^ appearing 
in the integral are even and periodic functions of 6. We 
recall that 

R = (p^+p'^ - 2ppcos 6 + (z ' - z) ^) 



and 



X = t - R/C 
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so that R, J'(p',z,t) and J'(p',z',x) are also even and 
periodic functions of 6. Finally U is independent of 3. 

We will concentrate on the integral over 3 in order to 
reduce it to its simplest expression. Since 

(j)' = (j) + 3 



we can write 



sin (j) ' = sin (p cos 3 + sin 3 cos <j) 

cos (|) ' = cos (j) cos 3 - sin (p sin 3 



and substitute this into the integral over 3 in 
equation (Cl). The integral over 3 for a complete period, of 
all terms that contain sin 3 will be zero since those terms 
are periodic and odd functions of 3. All the other terms will 
produce a non-zero integral and will be retained. 

After this simplification, equation (Cl) will become: 



i 1 ® 

sin (j) = -2H^ sin (p + sin (p / p'ds 



2tt , , ^ 



<(■ 



which after repressing the common term sin (p , 
final equation: 

J. = A'as’ 



0 R 



[J'F, cos 3 + VJ! sin^3] d3 . 

SI (p 



^3]d3 (C2) 

leads to the 



(C3) 
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Following the same reasoning for equation (3.5), and repressing 
the coininon term cos 4> leads to the following final equation: 



• 1 s 2 tt t "I T 

= -2H^ sin a + 5 ^ J p'ds' J 



. 2 , 



[J' Usings +J!F-cosS]dS (C4) 

S (p 2 

In both equations (C4) and (C3) the currents are: 



J_ A i(Pf2,t) 

s,9 s,4> 



s , <p s ,4) 



T being the retardated time t - R/C, and the incident field 
is : 



= H^(z,t) 
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OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOriOOOOOO 



SYS IN DD DATA 



APPENDIX D 



^ TRANSIENT RESPONSE OF A STRAIGHT THIN WIRE 

4c 

4e:jc:<e:<e4c4c4c4c4e4c4c4c4c4e4c4c4c4:4c3^4c:$:4e4c4e4c4e4c4c4c4:4c4c4c4c4e4c4c4c4c4c4c4c4c4c4c4c4c4c4c 



PURPOSE: 

3jc4c4c4c4c4c4c :jc 

THIS PROGRAM COMPUTES THE TRANSIENT TIME DOMAIN 
RESPONSE OF A STRAIGHT THIN WIRE. 

THE PROGRAM COMPUTES THE CURRENTS INDUCED ON THE 
WIRE 6Y AN INCIDENT IMPULSE FIELD, AND THE SCATTERED 
FIELD, FOR ANY ANGLES OF INCIDENCE AND SCATTERING. 



USAGE: 

THE PROGRAM CAN BE USED AFTER PROPER DEFINITION 
OF THE INITIAL PARAMETERS OF THE WIRE, AND DEFINITION 
OF THE DESIRED ANGLES OF INCIDENCE AND REFRACTION. 

THE WIRE IS ALIGNED WITH THE Z-AXIS STARTING AT Z=0 
THE ANGLES ARE MEASURED WITH RESPECT TO THE Z-AXIS 



PARAMETER DEFINITION 

Jfc 4(4: 4c:Ce 4c :0e 



8ETA(400,40) :ARRAY DEFINING THE CURRENTS FOR UP 

TO 40 SEGMENTSAND 400 TIME INCREMENTS. 

GAMMA! 400,40) : ARRAY DEFINING THECHARGES FOR UP 

TO 40 SEGMENTS AND 400 TIME INCREMENTS. 

NS: NUMBER OF SEGMENTS IN WHICH THE WIRE IS DIVIDED 

SL: total length OF THE WIRE IN METERS. 

A: RADIUS OF THE WIRE IN METERS. 

its: number OF TIME INCREMENTS. 

ITS MUST BE A MULTIPLE OF 80 FOR PLOTTING 
CONVENIENCE 

TI: ANGLE OF INCIDENCE IN DEGREES. 

TS: ANGLE OF REFRACTION IN DEGREES. 

DZ: LENGTH OF A SEGMENT (L/NS) 

DT: TIME INCREMENT IN SECONDS (DT=DZ/C) 

AA: GAUSSIAN IMPULSE PARAMETER. IN SEC-1. 
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TMAX: TIME DELAY OF THE IMPULSE FROM THE FIRST TIME 
OF INCIDENCE TO ITS PEAK. ( TYPICALLY 0T*NS/2» 

El; INCIDENT FIELD (GAUSSIAN PULSEK 

ES(400): ARRAY CONTAINING THE COMPUTED SCATTERED 
FIELD UP TO 400 TIME INCREMENTS. 

DELAZ : TIME DERIVATIVE OF THE MAGNETIC POTENTIAL. 

DELPHI; TIME DERIVATIVE OF THE ELECTRIC POTENTIAL. 

C; SPEED OF LIGHT. 

Cl; CONSTANT (MU/A^^'PIJ. 

C2; CONSTANT ( 1 /4*PI *EPS » . 

X(30»,Y(30»;C00RDINATES REQUIRED FOR SUBROUTINE 
» PLOTP* . 



REMARKS; 

THE ANGLES OF INCIDENCE AND REFRACTION DESIRED ARE 
ENTERED IN DATA CAROS ACCORDING TO FORMAT 2F6.1. 

EACH DATA CARO CONTAIN ONE ANGLE OF INCIDENCE 
AND ONE ANGLE OF REFRACTION. 

THE COMPUTATION IS DONE FOR AS MANY DATA INPUTS AS 
THERE ARE DATA CAROS. THE LAST DATA CARO MUST CONTAIN 
THE DATA NUMBER 1.0 IN ORDER TO TERMINATE THE PROGRAM. 



SUBROUTINES USED: 

3fl(3ic9fie:2&c TtCTbi sdcsdc :«dcsk :de^;6c :6c :dc 4c 

THE PROGRAM USES A LIBRARY SUBROUTINE NAMED PLOTP 
FOR THE PLOT OF CURRENTS AND FIELD AS A FUNCTION 
OF TIME. 



METHOD OF SOLUTION: 

THIS IS THE SOLUTION OF THE E.F.I.E. EQUATION 
APPLIED TO A THIN WIRE t ACCORD ING TO THE METHOD 
DESCRIBED BY HARRINGTON AND SAYRE IN 'TIME DOMAIN 
RADIATION AND SCATTERING BY THIN WIRES * (IN *APPL. 
SCI. RES. 26- SEPT 1972) 
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c 

c 

C J MAIM PROGRAM * 

C * * 

(^ ««4c4c4c4c4c4c4c4c4c4c4c4c 4( 4c4c 

c 

c 

c 

DIMENSION BETA(400,40I , GAMMA (400 »40), ESI 400) 

DIMENSION X(400) f Y(400) 

C 

C 

C CONSTANTS INITIALIZATION 

C 

C=3.0E+8 
Cl=l.0E-7 
C2=9.0E+9 
PI=3. 141592654 
C 
C 

C PARAMETERS INITIALIZATION 

C 

C COMPUTATION PARAMETERS 

C 

NS=40 

ITS=400 

NSl=NS-l 

C 

C WIRE PARAMETERS 

C 

SL=1.04 
A=7.75 E-4 
DZ=SL/NS 
OT=DZ/C 
C 

c 

C GAUSSIAN IMPULSE PARAMETERS 

C 

AA=0.67 S+9 
TMAX=DT*18.0 
C 

10 CONTINUE 
C 

C READING OF DESIRED ANGLE OF INCIDENCE TI 

C AND ANGLE OF REFRACTION TS. 

C 

READ(5,9510) TI,TS 
9510 FORMAT! 2F6. 1) 

IF(TI.EQ.l.O) GO TO 999 
C 

C WRITE HEADING 

C 

WRITE(6,960l) TI,TS 

9601 FORMAT! IHltlOX, ‘CURRENTS INDUCED WHEN FIELD INCIDENT*, 
2*AT ANGLE* ,IF6.1,//,10X,»FAR FIELD SCATTERED TOWARDS*, 
3 ‘ANGLE* , IF6. I ,//// , 8X, * TIME* , 15X , * CURRENT ( AMP) * ,/// I 
C 

C TRANSLATION FROM DEGREES TO RADIANS. 

C 

TIR=PI#TI/180.0 

TSR=PI^TS/180.0 

AI=GOS(TIR> 

AR=COS(TSR) 

C 

c 
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C START OF CURRENTS COMPUTATION 

C 

c 

Sl»0.5 
S2=-0.5 
S3=( A/0Z>**2 

FO=ALOG(C Sl + SQRTC Sl**2+ S3M /( S2-J-SQRT1 S2=«‘=«‘2+S3 ) I) 

C 

C INITIAL CONDITIONS SETTING 

C 

DO 100 L=lfNSl 
Ll=L-l 

IF(TI .LE. 90.01 GO TO 120 
L2=L-NS 
GO TO 130 
120 CONTINUE 
L2=L 

130 CONTINUE 

EI= EXP(( -AA*«2)*(DT*( 1-L2*AI )-TMAX)=<'=<‘2) *SIN (TIR) *3.0 
8ETA(1,U=EI/(2*C1<'F0/DT) 

IF(L.EQ.l) GO TO 110 

GAMMA <i,L )=-C 1.0/C )#( 8ETA( 1, U -SETA( 1,L1) J 
GO TO 100 
C 

C 8CUN0ARY CONDITIONS SETTING 

C 

110 CONTINUE 

GAMMA(1,1 )*-BETA (1»1 )/C 

100 CONTINUE 

GAMMA(1,NS)= BETA CltNSU/C 
C 

C START ALL CURRENTS COMPUTATION 

C AT EACH TIME INTERVALfANO FOR ALL SEGMENTS 

C 

DO 200 J=2fITS 
C 

DO 300 L*ltNSl 
C 

DELPHI =0.0 

OELAZ =0.0 

DO 310 K»1,NS1 
K1=IA8S(L-K) 

K2=J-K1 
SUKl+0.5 
S2=Kl-0.5 
S3=(A/DZ) **2 

F =ALOG(( Sl-i-SQRT( Sl^^Z+SS))/! S2<-SQRT( S2**2<-S3 ) ) ) 





IF(Kl.EQ.O) GO TO 


320 




IF(K2.LE.O) GO TO 


310 




OELAZ =DELAZ 


■••BETACK2,K)=(‘F 


320 


CONTINUE 

K3=K2-1 






IF(K3.LE.O) GO TO 


310 




OELAZ =DELAZ 


-BETA(K3,K)#F 


310 


CONTINUE 






OELAZ =DELAZ 


♦Cl 



DO 330 K=1,NS 
K2=IA8S CL-K» 

K1=IABS (L-K>H) 

K3=J-1-K1 
K4= J-1-K2 

IF(K3.LE.0» GO TO 340 

Sl=Kl+0.5 

S2=Kl-0.5 

F =ALOGC ( Sl+SQRTC Sl**2*S3)) /( S2+SQRT( S2*<'2 + S3 ) )) 
DELPHI ^DELPHI +GAMMA C K3 ♦ K) #F 

C 
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3^0 



330 



C 

c 

c 



c 

c 

c 

331 

332 



333 

334 



350 

300 



C 

C 

C 

c 

c 

9690 

200 

C 

C 

C 

C 

C 

C 

C 

c 



c 

c 

c 



c 

c 



CONTINUE 

IF(K4.LE.0» GO TO 330 

Sl = K2+0<,5 

S2=K2-0.5 

F =ALOGn SU-SQRT< Sl^*2<- S3 J ) /( S2+SQRT( S2*»2 + S3 )) ) 
DELPHI =0ELPHI -GAMMA(K4, K)*F 

CONTINUE 

DELPHI =OELPHI *C2 

IFITI.LE. 90.01 GO TO 331 

FOR ANGLES OF INCIDENCE MORE THAN 90 DEGREES 



L2=L-NS 
GO TO 332 

FOR ANGLES OF INCIDENCE LESS THAN 90 DEGREES 

CONTINUE 
L2 = L 

CONTINUE 

IF (J.GT.50 ) GO TO 333 

EI= £XP< (-AA#>!‘2)«=( J*0T*(l-L2=«'An-TMAX)**2» *SIN(TIRJ*3 
GO TO 334 
CONTINUE 
EI*0.0 
CONTINUE 

BETA<J,U*( EI-OELAZ 
J1=J-1 
L1=L-1 

IFIL.GT.l ) GO TO 350 
GAMrtA(J,U= GAMHA(J1,L 
GO TO 300 
CONTINUE 

GAHMA(JtL)= GAMMA! J1,L 
CONTINUE 

GAMMA! J, NS) = GAMMA ! Jl , NS ) +BET A ! J , NS 1 J/C 
KK=NS/2 



/OT-OELPHI 

)-BETA!J,L 

)-BETA!J,L 



/DZ)/!C1«F0/0T) 



)/C 



)/C +BETA! J,L1)/C 



WRITE THE CURRENT INDUCED AT THE CENTER OF THE WIRE 



WRITE!6,9690) J,BETA!J,KK) 
FORMAT!10X,I4f21X,E12.5) 

CONTINUE 

END OF CURRENTS COMPUTATION 



START OF SCATTERED FIELD COMPUTATION 



00 400 LP=lfITS 
ES!LP)=0.0 
DO 500 K=1,NS1 
IF !TS.GT.90.0) GO TO 501 
IF ITS. EQ. 90.0) GO TO 503 

FOR REFRACTION ANGLE LESS THAN 90 DEGREES 

KP^ K 

NQ=LP-1-! KP-0.5)*AR 
NQP=LP-1- !KP+0.5)*AR 
GO TO 502 
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C FOR REFRACTION ANGLE MORE THAN 90 DEGREES 

C 

501 CONTINUE 

kp=k-ns 

NQP=LP-1-(KP-0.5)«'AR 
NQ=LP-L-( KP-»-0.5)#AR 
GO TO 502 
C 

C FOR REFRACTION ANGLE EQ. 90 DEGREES. 

C 

503 CONTINUE 
NQ=LP-l 
NQP=NQ 

502 CONTINUE 
NQl=NQ-l 
NQ2= NQP-1 

IF(NQ.LE.O) GO TO 500 
IF(NQ.EQ.l) GO TO 510 
IF( NQ.EQ.NQP) GO TO 520 
ZQ=(LP-1-NQ)/AR •►0.5 
EPS=D2*(ZQ-KP) 

DEL=OZ-EPS 

ES(LP)=ES(LP) •►£?$*( BETA! NQfK) -BETA! NQ1,K) ) 

IFCNQ.EQ. 2) GO TO 530 

ES(LP)=ES(LP) ►DEL* (BETAINQP, K)-BETA( NQ2, K ) ) 

GO TO 500 
510 CONTINUE 

ES(LP)=ES(LP)+0Z*2* B£TA(NQ,K) 

GO TO 500 
520 CONTINUE 

ES(LP)=ES (LP)+OZ*(BETA(NQtKi-BETA(NQl,K) ) 

GO TO 500 
530 CONTINUE 

ES<LP)=ES(LP)+0Z*2* BETAINQP, K) 

500 CONTINUE 

CT=-SIN(TSR)*a.O E-7)/DT 
ES(LP)=CT#ES< LP) 

400 CONTINUE 
C 
C 
C 

C END OF FIELD COMPUTATION 

C 

C 

c 

C PLOT OF CURRENT AT CENTER OF THE WIRE 

C PLOT EVERY MM POINT. 

C 

MM=ITS/100 
00 600 I=MM,ITS,MM 
K=I/MM 
X(K) =I=«'OZ 
KK=NS/2 

Y(K)= BETAIIfKK) *1000.0 
600 CONTINUE 

MM= ITS/80 
WRITE(6,9630) 

9630 FORMAT! IHl) 

MO0CUR=0 

CALL PL0TPIX,Y,80 ,MODCUR» 

WRITE (6,9600) TI 

9600 FORMAT!// //,5X, 'CURRENT INDUCED AT CENTER OF WIRE',' 

2 SCA'^'TERER' , /,5X, 'FIELD INCIDENT AT ANGL E* , 1F6 . 1,/ // , 

3 lOX, 'X-AXIS TIME IN LIGHT-METERS', 

4 //,10X,' Y-AXIS CURRENT IN MILLIAMPS') 

C 
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c 

C PLOT OF SCATTERED FIELD 

C 

DO 700 I=MMtITS,MM 
K=I/MM 

Y(K)=£Sm ♦1000.0 
700 CONTINUE 

WRITE (6,9640) 

9640 FORMAT! IHl) 

M0DCUR=0 

CALL PLOTP(X,Y,80 ,M0DCUR) 

WRITE (6,9620) TS 

9620 FORMAT!////, 5X, 'NORMALIZED FAR FIELD BACK SCATTERED',' 
2BY WIRE', /,5X,'AT ANGLE' , 1F6 .1 , / / , lOX, ' X-AX I S TIME',' 
3 LIGHT METERS* ,// ,10X, 'Y-AXIS FAR FIELD MILLIVOLTS') 

C 

C 

C END OF PROGRAM GO BACK TO NEXT ANGLE IF ANY 

C 

GO TO 10 
999 CONTINUE 
STOP 
END 
C 

c 

c 

c 

c 

c 

C EXAMPLE OF DATA CARDS 

C 
C 
C 

90.0 90.0 

l.O 
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APPENDIX E 



:(c 4c 4c 4c 4c 4c 4c 4c 4c :(c 4c4c 4c 4c4c 4: C 4c 4c 4c 4c 4c 4c 4: 4c 4c 4c 4c 4c # 4c # 4c 4c 4c 4c 4c 



♦ TIME DOMAIN RESPONSE AND INVERSE SCATTERING J 

4c 4c 

=)« FOR A PERFECT CONDUCTOR 30DY OF REVOLUTION * 

4c 4c 



4c 4c 4c 4<4c 4c 4c 4c 4c 4c 4c 4c4c 4c4c 4c 4c 4c 4c 4c4c 4c 4c 4c4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4: 4c 4c 4c 4c 4: 4c 4c 4c 4c 4c 4c4c 4c 



NAME OF PROGRAM : GAIL 

4c 4c4c4c 4c 4c4c4c4c 4c 4c 4c4c 4c 4c 4c 



PURPOSE: 

4c 4c 4c 4c4c 4c 4c 4c 

THIS PROGRAM COMPUTES THE TRANSIENT TIME DOMAIN 
RESPONSE OF A PERFECT CONDUCTOR BODY OF REVOLUTION, 
AND APPLIES IT TO THE SOLUTION OF THE INVERSE 
SCATTERING PROBLEM, (CCMPUTATI ON OF THE SHAPE OF THE 
BODY FROM ITS TIME DOMAIN RAMP RESPONSE) 
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usage: 



THE PROGRAM CAM BE USED IN ONE OF FOUR WAYS 
DEPENDING OF THE DESIRED OUTPUT AND CONSEQUENTLY OF 
THE VALUE ASSIGNED TO THE INITIATING PARAMETERS 'INC* 
AND ’INVER*. 

THE POSSIBLE USES ARE: 

1. COMPUTATION OF THE TRANSIENT TIME DOMAIN IMPULSE 
RESPONSE (BACKSCATTEREO FIELD ) tWHEN THE INCIDENT 
FIELD IS A GAUSSIAN IMPULSE .(INC=1) 

THE CURRENTS INDUCED ON THE BODY CAN BE MADE AVAILABLE 
AT THE OUTPUT IF DESIRED (FOR EMP PROBLEMS). 

THE INPUTS REQUIRED ARE: 

A. THE SPACE SAMPLES OF THE CONTOUR FUNCTION OF 
THE BODY,RO(Z),ENTERD IN DATA CAROS. 

B. THE NUMBER OF SEGMENTS SUBDIVIDING THE CONTOUR 
OF THE BODY (NS<20). 

C. THE LENGTH OF A TIME INCREMENT (DT). 

0. THE NUMBER OF TIME SAMPLES NEEDED ( ITS<= 60). 

E. THE LENGTH OF EACH SEGMENT DS(20). 

F. THE PARAMETERS DEFINING THE GAUSSIAN IMPULSE. 

2. COMPUTATION OF THE RAMP RESPONSE OF THE BODY 
WHEN THE INCIDENT FIELD IS A RAMP.(INC=2 AND INVER= I ) 

THE INPUTS REQUIRED ARE: 

SAME AS A THRU E LISTED ABOVE. 

3. INVERSE SCATTERING SOLUTION STARTING FROM THE 
BODY ITSELF. (THIS CASE IS A SIMULATION PROCEDURE, 
WHERE THE PROGRAM RETRIEVES BY INVERSE SCATTERING 
THE SHAPE OF THE KNOWN BODY. (INC=2 AND INVER=l). 

THE INPUTS REQUIRED ARE: 

SAME AS A THRU E LISTED ABOVE, PLUS THE LENGTH OF 
A SEGMENT FOR NEW SUBDIVISION OF THE BODY.(OSO) 

4. INVERSE SCATTERING SOLUTION STARTING FROM THE 
MEASURED GIVEN RAMP RESPNSE OF AN UNKNOWN BODY (HSFMI. 
THIS CASE IS THE TRUE COMPUTATION OF THE SHAPE OF THE 
BODY WHEN ONLY ITS RAMP RESPONSE IS KNOWN. 

(INC=2 AND INVER=2I 

THE INPUTS REQUIRED ARE: 

A. THE TIME SAMPLES OF THE RAMP RESPONSE . (ENTERED 
IN DATA CAROS) 

8. THE LENGTH OF THE DESIRED SEGMENTS. (DSO) 

C. B THRU 0 LISTED ABOVE. 



NOTE: NO INVERSE SCATTERING IS COMPUTED WHEN INC=l . 
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PARAMETER DEFINITION: 



INC: DEFINES THE TYPE OF INCIDENT FIELD 

INC=1 THE INCIDENT FIELD IS A GAUSSIAN IMPULSE. 
INC=2 THE INCIDENT FIELD IS A RAMP. 



IHVER : DEFINES THE TYPE OF PROGRAM REQUIRED 

INVER=0 NO INVERSE SCATTERING 

INVER=l THE PROGRAM COMPUTES THE INVERSE SCATTERING, 
STARTING FROM A GIVEN BODY SHAPE. 

INVER=2 THE PROGRAM COMPUTES THE INVERSE SCATTERING 
FROM A GIVEN RAMP RESPONSE OF THE BODY. 



NS: NUMBER OF SEGMENTS SUBDIVIDING THE CONTOUR (<=20) 

ITS: NUMBER OF TIME SAMPLES NEEDED ( <= 60) 

DT: LENGTH OF THE TIME INCREMENT (SECONDS) 

DS(20) : LENGTH OF THE SEGMENTS SUBDIVIDING THE 

CONTOUR OF THE BODY. ( MAXIMUM 20 SEGMENTS) 

TO BE ENTERD IN DATA STATEMENT. 

Z(21): Z-COORDINATE OF THE SPACE SAMPLE POINTS ON THE 
CONTOUR OF THE BODY. ( MAXIMUM 21 POINTS ENTERED 
DELIMITING 20 SEGMENTS ).( METER) 

RO (21) :RADIUS of THE BODY AT THE SPACE SAMPLE POINTS 
DEFINED BY Z.(MAX. 21 ). (METER) 

CJS (20,60): SPACE TIME SAMPLES OF THE S COMPONENT OF 
THE CURRENTS INDUCED ON THE BODY. (UP TO 
20 SPACE SAMPLES AND 60 TIME SAMPLES EACH) 

CJPHK20, 60) : SPACE TIME SAMPLES OF THE PHI COMPONENT 
OF THE CURRENTS INDUCED ON THE BODY. (UP TO 
20 SPACE SAMPLES AND 60 TIME SAMPLES EACH) 

HSF (60): TIME SAMPLES OF THE SCATTERED FIELD. UP TO 

60 SAMPLES. 

HSFM(60): TIME SAMPLES OF THE INITIAL RAMP RESPONSE. 

( THE MEASURED RAMP RESPONSE) 

HSFC(60): TIME SAMPLES OF THE COMPUTED RAMP RESPONSE 

OF THE BODIES OBTAINED BY INVERSE SCATTERING 

NP(20): NUMBER OF PATCHES ON THE RINGS DEFINED BY 
EACH SEGMENT. 

B(20): ANGLE SUSTAINED BY THE PATCHES ON EACH RING. 

ALFAK: SLOPE OF THE CONTOUR OF THE BODY CONTOUR 
FUNCTION , MEASURED AT THE CENTER OF THE 
EACH SEGMENT. 

SI (20): S1N( ALFAK) FOR EACH SEGMENT. 

CO (20): COS( ALFAK) FOR EACH SEGMENT. 

EPS: NORMALIZED MEAN SQUARED ERROR ON THE BODY SHAPE 
COMPUTED BY INVERSE SCATTERING PROCEDURE. 
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ITER: NUMBER OF ITERATIONS EXECUTED FOR INVERSE 
SCATTERING SOLUTION. 

C: SPEED OF LIGHT (METER/SEC) 

CP : CONSTANT C^'OT. 

OTHER PARAMETERS USED ARE DEFINED IN THE PROPER 
SUBROUTINE. 



REMARKS : 

1. THE SHAPE OF THE BODY IS ENTERD IN DATA 
CAROS WITH THE FORMAT 2F10. 5 

2. THE TIME SAMPLES OF THE MEASURED RAMP 
RESPONSE HSFM ARE ENTERD IN DATA CAROS WITH THE 
FORMAT 1F10.5 

THE VALUES OF THE SAMPLES MUST BE NORMALIZED, THAT 
MEANS MULTIPLIED BY C AND R (RANGE TO TARGET) 

3. IN THIS PROGRAM THE INCIDENT FIELD H IS 
ASSUMED TO BE VERTICAL (HORIZONTAL POLARIZATION). 
THE DIFFERENCE BETWEEN THE RESPONSE TO THE TWO 
POLARIZATIONS IS ONLY THE SIGN OF THE RESPONSE. 

( SEE ALSO REMARKS IN SUBROUTINE SHAPE) 



SUBROUTINES USED: 

1. SUBROUTINE 'FIELD* : COMPUTES THE BACK SCATTERED 
FIELD FOR BOTH RAMP AND IMPULSE EXCITATION. 

2. SUBROUTINE 'SHAPE* : COMPUTES WITH A METHOD OF 
ITERATIONS. THE SHAPE OF THE BODY FROM ITS RAMP 
RESPONSE. (SEE SUBROUTINE SHAPE FOR DETAILS.) 



3. LIBRARY SUBROUTINE 'PLOTP* FOR PLOTS OF FIELDS 
AND SHAPE OF THE BODY. 



METHOD OF SOLUTION: 

<c «4c4c :;c 4a(c 4cc»: « « 4c« 4c« 4c 4c 4c4c 

SEE SUBROUTINE FIELD FOR THE COMPUTATION OF IMPULSE 
AND RAMP RESPONSES. 

SEE SUBROUTINE SHAPE FOR THE COMPUTATION OF THE 
SHAPE OF THE BODY. 

INVERSE SCATTERING: THE FIRST APPROXIMATION OF THE 
SHAPE IS OBTAINED FROM ITS MEASURED RAMP RESPONSE. 

A NEW RAMP RESPONSE IS COMPUTED FOR THAT APPROXIMATE 
BODY. THE PHYSICAL OPTICS PART OF THE RESPONSE IS 
EXTRACTED, AND IS USED TO COMPUTE A NEW SHAPE, AND SO 
ON , UNTILL A MINIMUM MEAN SQUARED ERROR IS OBTAINED 
BETWEEN THE LAST COMPUTED RESPONSE AND THE MEASURED. 
(SEE ALSO SECTION 5 OF THIS THESIS) 
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Q a}c:^c:Jr:4c:Cc:^K:!#:3t:3jc3}c:?:Cc 5js:{c <s :? ajc 

C ^ ^ 

C MAIN PROGRAM =<' 

C * * 

C 

C SPHERE EXAMPLE 

C 

DIMENSION X(60) 

DIMENSION CJS(20t60 ) ,CJPHI ( 20 ,60) » HSF(60) 

DIMENSION Z(21) ,R0( 21) ,NP(20) 

DIMENSION SI (20) , CO ( 20 ) , B (20) ,DS (20) 
DIMENSION HSFM(60) ,HSFC (60) 

COMMON HSF,Z,RO,ITS ,NS,NS1 

COMMON XMAX 

DATA HSFM/60*0.0/ 

DATA HSFC/60*0.0/ 

DATA DS/20=«'0. 20900/ 

C 

C 

C START PROGRAM 

C 

C DEFINE TYPE OF PROGRAM, INVER AND INC 

C 

INC=2 

INVER=1 

C 

C COMPUTATION PARAMETERS 

C 

NS = 15 
NS1=NS+1 
ITS=60 
IT=ITS 
DSO=0.209 
C=3.0E+8 
0T=0.2/C 
CP=C*DT 
C 

IF (INC.EQ.l) GO TO 11 
IF (INVER .EQ.2) GO TO 15 
C 

C STARTING POINT FOR INV£R=1 AND INVER=0 

C 

11 CONTINUE 

C 

C READ THE SHAPE OF THE BODY 

C 

DO 10 K=1,NS1 
READ(5,9510) Z(K),RO(K) 

9510 F0RMAT(2F10.5) 

10 CONTINUE 

C 

C PLOT THE INITIAL SHAPE OF THE BODY 

C 

NS2=NS1+1 

XMAX = Z (NS1)*1.25 
Z(NS2)= XMAX 
RO (NS2)= Z(NS2) 

WRITE(6,9620) 

9620 FORMAT UHl) 

MODCUR=0 

CALL PLOTP( Z,R0,NS2,M00CUR) 

WRITE (6,9621) 

9621 FORMAT ( /// , 20X , » I NIT I AL SHAPE OF BODY* , // ,2 IX , 
2 'X-AXIS Z, Y-AXIS RADIUS (RO)') 

C 

C 
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C COMPUTE THE BACK- SCATTERED FIELD 

C 

C 

CALL FIELD (OS,CP , OT , INC I 
C 
C 

C END OF PROGRAM FOR INVER=0 AND INVER=l 

C 

IF( INC.EQ.l) GO TO 99 
IF (INVER .EQ-O) GO TO 99 
IF (INVER .EQ.U GO TO 19 
15 CONTINUE 

C 

C STARTING POINT FOR INVER=2 

C 

INC=2 

C 

C READ THE GIVEN RAMP RESPONSE 

C 

DO 12 N=l,ITS 

READ (5,9520) HSFM (N» 

9520 FORMAT (1F10.5) 

HSF(N)=HSFM(N) 

12 CONTINUE 

GO TO 920 
C 

19 CONTINUE 

C 

C STORE THE COMPUTED RAMP RESPONSE 

C 

DO 910 N*1,ITS 
HSFM(N)=HSF(N) 

910 CONTINUE 
920 CONTINUE 
C 

C STARTING POINT FOR INVERSE SCATTERING COMPUTATION 

ITER=1 

C 

C COMPUTATION OF THE FIRST APPROXIMATION OF THE SHAPE 

C 
C 
C 

CALL SHAPE ( OS , ITER,CP ,DSO , I T) 

C 

C 

C 

C COMPUTE SUM OF SQUARES OF THE TIME SAMPLES OF THE 

C ORIGINAL RAMP RESPONSE. 

C 

ITS=IT+10 
SIGMA1=0. 0 
DO 940 N=l, ITS 
SIGMA1=SIGMA1+HSFM(N)**2 
940 CONTINUE 

C 

950 CONTINUE 

C 

C COMPUTE THE RESPONSE OF THE APPROXIMATE BODY. 

C 

C 

C 

CALL FIELD ( DS,CP , DT , INC ) 

C 

C 

C 

C 
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c 

c 



960 

C 

C 

c 

9666 



961 

C 

C 

C 



C 

C 

C 

c 

c 

c 

9649 



9651 

970 

C 

c 

c 

c 

c 

r 

c 

c 

999 

C 

c 

c 

c 

9625 

99 

C 

C 

c 



c 



COMPUTE NORMALIZED MEAN SQUARED ERROR. 

SIGMA2=0. 0 
00 960 N=1,ITS 

SIGMA2=SI GMA2+(HSFM(NI-HSF(N) >**2 
CONTINUE 

EPS=SIGMA2/SIGMA1 
WRITE ERROR 
WRITE (6,9666) EPS 

FORMAT i //,15X, 'NORMALIZED MEAN SQ. ERROR » ,1F10 . 5 1 
IF (ITER.GT.l) GO TO 961 
EPS1=1. 

CONTINUE 

END THE PROGRAM WHEN THE ERROR IS MINIMUM 

IF(EPS.GT.EPSl) GO TO 999 
EPS1=EPS 

CONTINUE WHEN THE ERROR IS NOT YET MINIMIZED 
ITER=ITER+1 

EXTRACTION OF THE PHYSICAL OPTICS RESPONSE 
WRITE (6,9649) 

FORMAT (IHl, lOX, 'PHYSICAL OPTICS RESPONSE ',//, 9X ,' N' , 
2 10X,'H( N) ' ,//) 

DO 970 N=1,ITS 
HSFC(N)=HSFC(N)+HSF(N) 

HSF(N)=ITER*HSFM(N)-HSFC(N) 

WRITE (6,9651) N,HSF(N) 

F0RMAT(5X,I5,5X,E12.5) 

CONTINUE 

USE PHYSICAL OPTICS RESPONSE TO COMPUTE A NEW SHAPE 



CALL SHAPE ( DS , I TER,CP ,DSO, I T ) 



RETURN TO COMPUTE THE RESPONSE TO THE NEW BODY 

GO TO 950 
CONTINUE 

END OF PROGRAM WHEN MINIMUM ERROR ACQUIRED 
WRITE FINAL ERROR 

WRITE (6,9625) ITER.EPSl 

FORMAT (//,20X, 'FINAL SHAPE AFTER ', 14, ' ITERATIONS ' ,/// 
2 ,20X, 'NORMALIZED MEAN SQUARED ERROR IS',1F10.5) 

CONTINUE 

END OF PROGRAM 

STOP 

END 
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:0c :flc :jc4e Jjc:iK:jc:je 4c 

4c 4e 

« SUBROUTINE FIELD ♦ 

4c 4 k 

4c4c 4c4c 4 kj4k 4k 30c4t 4:4c 4c4c :^r4K4c3iK 4c 4K30C4C 30: 4e 



PURPOSE: 



THIS SUBROUTINE PROVIDES THE SOLUTION TO THE 
MAGNETIC FIELD INTEGRAL EQUATION APPLIED TO A 
BODY OF REVOLUTIONfEXCITED BY A GAUSSIAN PULSE 
OR A RAMP FIELDS. THE CURRENTS INDUCED ON THE BODY 
AND THE BACK-SCATTERED FIELD ARE COMPUTED. 

usage: 



CALL FIELD (OS , CP t DT, INCI 



PARAMETER DEFINITION: 
4c :jc « « 4: «4 e « 4= :4c 4c 3ic^:{c<c 3jc 4c 



HI: INCIDENT FIELD 
C3»C4: CONSTANTS. 



METHOD OF SOLUTION: 

4: 4> 4c4c 4c 4c4c 4c4c 4c 4c 4c4c 4c4c 4c 4c 4c 4c 



THE METHOD OF SOLUTION IS DESCRIBED IN SECTION 4 
OF THIS THESIS. 
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c 

c 



SUBROUTINE FIELD I D S, CP ,DT, INC J 



DIMENSION X(60) 

DIMENSION CJS(20,60 ) , CJ PHI ( 20 » 60 } , HSF(60) 

DIMENSION Z(2L) ,R0( 21) ,NP(20> 

DIMENSION SI (20) , CO ( 20 ) , B (20) ,DS (20) 

COMMON HSF,Z,R0, ITS,NS,NSl 

CONSTANT DEFINITION 

PI=3. 141592654 
C=3.0E+8 

C4=-1./(DT*4.0=«'C) 

C3=l./( PI*2.0) 

GAUSSIAN IMPULSE PARAMETERS 

AA=0.6 E+9 
TMAX=DT^5.0 
A2=AA*AA 

00 20 K=1,NS 
Ki=K+l 

COMPUTE THE SLOPE OF THE SEGMENTS 

WK=(R0(K1 )-R0(K))/(Z(Kl)-Z(K) ) 

ALFAK=ATAN( WK) 

SI( K)=SIN(ALFAK) 

CC( K)=C0S(ALFAK) 

COMPUTE THE CENTRAL SPACE SAMPLE POINTS 

Z(K)=(Z(K1)+Z(K) )/2.0 

R0(K)=( R0(Kl)-i-R0( K) )/2.0 

COMPUTE THE NUMBER OF PATCHES IN EACH RING 

PK= PI4=R0( K)/DS( K) 

NP(K)=2=«'IFIX( PK) 

IF (NP(K) .GT.O) GO TO 23 
NP( K)=l 

COMPUTE THE ANGLE SUSTAINED BY THE PATCHES OF EACH 
RING. 

CONTINUE 

B(K)=2.0*PI/NP(K) 

INITIAL CONDITIONS OF THE CURRENTS 

TZK=1.0+Z(K)/CP 
DO 25 N=l,2 

IF (N.LT.TZK) GO TO 26 
IF( INC.EQ.l ) GO TO 24 

WHEN INCIDENT FIELD IS A RAMP 



HI=(N-1) *CP-Z(K) 

GO TO 27 
4 CONTINUE 

WHEN INCIDENT FIELD IS AN IMPULSE 

ARG=(N-1) *DT-TMAX-Z(K)/C 
HI = (N-l)=<fEXP(-A2*ARG=J'ARG) 
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26 

27 

C 

C 

c 



25 

20 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 



c 



410 

420 

C 

C 

c 



c 

c 

c 

c 

c 

c 

c 



c 

c 



GO TO 27 
CONTINUE 
HI=0.0 
CONTINUE 

INITIAL CURRENTS INDUCED 

CJPHl(KtN)=-2.0^HI*SI(K) 

CJS(KtN)=-2.0*HI 

CONTINUE 

CONTINUE 



COMPUTE THE CURRENT COMPONENTS 

COMPUTE FIRST THE INTEGRAL PART OF THE EQUATION 
AT TIME N 

00 100 N=3fITS 

AT THE K-TH SAMPLE 

DO 200 K=1,NS 

CJS(K,N)=0.0 

CJPHI(K,N)=0.0 

INFLUENCE OF ALL OTHER SPACE CURRENTS 

DO 300 1=1, NS 
NP2=NPU) 

DO 400 M=1,NP2 
IF (M.EQ.l) GO TO 410 
GO TO 420 

CONTINUE 

IF ( I .EQ.K) GO TO 400 
CONTINUE 

COMPUTE DISTANCE BETWEEN SAMPLES 

CM=COS( (M-l)*Bn» ) 

SM=SIN( (M-IM'B(I) ) 

V=Z(I)-Z( K) 

R2= ( ROC K) <=R0( K) ♦ROC I ) *R0 C I ) +V*V-2. 0*R0C I » ♦RO ( K) =fl‘CM) 
R=SQRT(R2) 

COMPUTE RETARDATION TIME 
Q=N-R/CP 

INTERPOLATION OF THE CURRENT AND THEIR TIME 
DERIVATIVE AT THE RETARDED TIME 

NQ=IFIX(Q) 

NQP=NQ+1 

NQM=NQ-1 

NQ2*NQ+2 

T1=Q-NQ2 

T2=Q-NQP 

T3=Q-NQ 

T4=Q-NQM 

IF CQ,LE.1.0) GO TO 400 

IF CQ.LE.2.0) GO TO 435 

IF (N.EQ.NQP) GO TO 430 

IF (N.EQ.NQ2) GO TO 432 
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FCUR POINT LAGRANGIAN INTERPOLATION 



C 

C 

C 



A30 

C 

c 

c 

c 

c 



A32 

C 

C 

c 

c 

c 



A35 



C 

c 

c 



A3 3 
C 
C 
C 



c 

440 



yi=T2#T3#T4*CJS( ItNQ2) /6.0-T3=«'T4*Tl*CJSI I ,NQP 1/2.0 
2 +T4*T1»T2*CJS( I » NQI /2.0-Tl*T2*T3*CJS( I »NOM)/6.0 

Y2=T2*T3*T4*CJPHI( I,NQ2)/6.0-T3*T4*Tl=«'CJPHI( I,NQP)/2.0 
2 +T4*TI«T2=«'CJPHI ( I,NQ 1/2. 0-Tl*T2*T3=!'CJPHI ( I *NQM)/6.0 
DJS=(T2=«'T3+T2*T4<-T3*T4I=«‘CJS( I,NQ2)/6*0 

2 -(T3#T4+T4<‘TI+T1»T3)*CJS( I,NQP)/2.0 

3 +{T4=«'TH-T4=!'T2+T1=<'T2)*CJS( I»NQ) /2.0 

4 -(Tl*T3 + T2=<‘T3<-T2*Tl)=«'CJS( I,NQM)/6.0 
DJPHI = (T2*T3+T2>^T4+T3*T4)*CJPHI( I ,NQ2)/6.0 

2 -(T3#T4<-T4*T1<-T1*T3)*CJPHIC I ,NQP)/2.0 

3 +(T4*TI+T4*T2 + T1=«‘T2)^CJPHI{ I,NQ)/2.0 

4 -(T1#T3+T2=«‘T3+T2=<'T1I*CJPHI(I ,NQM»/6.0 
GO TO 440 

CONTINUE 



TWO POINT LAGRANGIAN INTERPOLATION 



Y1=T4»CJSU,NQ l-TS^CJSdfNQM J 

_ 



ti^jy )-T3«CJPHI (I ,NQM) 
)-CJS( I,NQMI 



Y2=T4*CJPHI(I 
DJS=CJS(IfNQ 
DJPHI=CJPHI(I ,NQ )-CJPHIU,NQM) 
GO TO 440 
CONTINUE 



THREE POINT LAGRANGIAN INTERPOLATION 



Y1 = T4*T3#CJS( I,'4QP)/2.0-T2*T4#CJS(ItNQ ) 

2 ■►T2*T3*CJS{ ItNQM) /2.0 

Y2=T4*T3#CJPHI( I,NQP)/2.0-T2=«'T4*CJPHI{I,NQ I 
2 +T2*T3*CJPHl(I»NQM)/2.0 

0JS = (T4-»-T3)^CJS(I,NQPJ/2.0-{T2+T4l*CJS( ItNQI 
2 +{T2*T3)=«'CJS{ l,NQM)/2.0 

DJPHI=(T4<-T3) ♦CJPHI ( I , NQP > /2. 0- ( T2+T4) *C JPHI ( I,NQ I 
2 +(T2+T3)=«‘CJPHI(I,NQM)/2.0 

GO TO 440 
CONTINUE 

IF (N.EQ.NQP) GO TO 400 
IF (N.EQ.NQ2) GO TO 433 

THREE POINT LAGRANGIAN INTERPOLATION 

Yl = T2*T3*CJS( I,NQ2)/2.0-T3=«'Tl=!'CJS{ ItNQP) 

2 +T2=!'Tl»CJSU,NQ)/2.0 

Y2=T2*T3*CJPHI( I,NQ2)/2.0-T3=>Tl=«'CJPHI( I,NQP) 

2 +T2=«'T1*CJPHI( I,NQ)/2.0 

0JS=(T2<-T3)*CJS( I,NQ2)/2.0-(T3<-Tl )*CJS< I,NQP) 

2 +(T2<-Tl )=<‘CJS( I,NQ)/2.0 

0JPHI=(T2+T3) *CJPHI(I,NQ2)/2.0-(T3+TlH'CJPHI{ I,NQP) 
2 •►(T2-«-Tl)*CJPHI{I,NQ)/2.0 

GO TO 440 
CONTINUE 

TWO POINT LAGRANGIAN INTERPOLATION 

Yl = T3=«'CJS( ItNQP)-T2*CJS(I,NQ) 

Y2=T3*CJPHI ( I ,NQP )-T2=«‘CJPHI ( I ,NQ) 

DJS=CJS(I »NQP )-CJS( I,NQ) 

DJPHI=CJPHI (I ,NQP )-CJPHI( I,NQ) 

CONTINUE 
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c 

C FINAL COMPUTATION OF THE INTEGRAL 

C 

Fl = V*CM=«'SI( I )+(RO(KI-RO(I l*CMI*CO( I I 

Xl = Fl=i‘CM 

X2=V=('(SM*SMI 

U=V*SI( D^SKKI+ROIKI^COIKI^SK II-RO( n*SI(K»<'CO( II 
F2=(R0( K)*CM-RO( 1 1 )=«'CO( KH-CM^C'V^S I ( Kl 
Yll=Yl*Xl 
Y22=Y2*X2 
XPl=F2*CM 
XP2=U*( SM*SM I 
YPI=Y2#XP1 
YP2=Yl#XP2 

Cl=OS( n»RO( II*B( I l/tCP*R2l 
R3=R*R2 

C2=0S( II* RO( I l*B(II/R3 
CJPHI(KfN)=CJPHI<K,NI+Cl*(DJPHI*XPl+0JS*XP2l 
2 +C2*( YPH-YP2I 

CJS(K,NI= CJS(K,NH-C1*(9JPHI*X2 >OJS*Xl H-C2* ( Y11+Y22 1 
400 CONTINUE 

C 

300 CONTINUE 

CJS(K,Nl=CJS(K»NI*C3 
CJPHK KfNI=CJPHI (KfNI *C3 
C 
C 

C END OF COMPUTATION OF THE INTEGRAL 

C 

C 

C ADDITION OF THE CURRENT INDUCED BY THE INCIDENT 

C FIELD AT THE COMPUTATION PO INT .( PHYSICAL OPTICS PART 

C 

TZK=1.0+2(KI/CP 
TD=TZK>2.*TMAX/DT 
IF (N.LT.TZKI GO TO 200 
IF( INC.EQ.ll GO TO 240 
C 

C WHEN INCIDENT FIELD IS A RAMP 

C 

HI=(N-ll*CP-Z(Kl 
GO TO 250 
240 CONTINUE 

IF (N.GT.TOI GO TO 200 
C 

C WHEN INCIDENT FIELD IS AN IMPULSE 

C 

ARG=(N-1I *OT-TMAX-Z(KI/C 
HI=EXP(-A2*ARG*ARGI 
250 CONTINUE 
C 

C FINAL CURRENT AT THE K-TH SPACE SAMPLE AT TIME N 

C 

CJS(K,NI=CJS(Kf NI-2.0*HI 

CJPHK X,NI=CJPHI(K,N l-2.0*HI#SI(KI 

C 

C NEXT SPACE SAMPLE 

C 

200 CONTINUE 

C 

C NEXT TIME 

C 

100 CONTINUE 
C 

C END OF COMPUTATION OF ALL SPACE-TIME CURRENTS 

C 
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IF UNC.eQ.ll GO TO 110 
WRITE(6t9623 ) 

9623 FORMAT! IHl ,10X, 'NORMALIZED FIELD (H*C*R)'f//f 
2 lOX, 'RAMP RESPONSE* ,//#9X# 'N* , 10X,'H(NM ♦//) 

GO TO 111 
110 CONTINUE 

WRITE(6#9723) 

9723 FORMAT! IHlflOXf'NORMALIZED FIELD (H*R) 

2 lOX, » IMPULSE RESPONSE' ,//,9X,'N' ,10X,'H!N) • ,//) 



START COMPUTATION OF THE BACK- SCATTERED FIELD 



I I CONTINUE 

AT THE N-TH TIME 

DO 500 N=l,ITS 
HSF!N)=0.0 

IF !N.EQ. 1) GO TO 650 
FROM EACH SAMPLE 
DO 600 K=1,NS 

COMPUTATION OF THE RETARDATION TIME 

g=N-Z(K)/CP 

NQ=IFIX(Q) 

INTERPOLATION OF THE TIME DERIVATIVE OF THE CURRENT 
AT THE RETARDED TIME 

NQP=NQ-H 

NQM=NQ-1 

NQ2=NQ+2 

T1=Q-NQ2 

T2=G-NQP 

T3=Q-NQ 

T4=Q-NQM 

IF !Q.LE. l.O) GO TO 600 
IF <Q.LE.2.0) GO TO 601 

IF !N.EQ.NQP) GO TO 602 
IF !N.EQ.NQ2) GO TO 603 



FOUR POINT LAGRANGIAN INTERPOLATION 

DJS=(T2<'T3+T2=*T4+T3*T4I *C JS ( K #NQ2 > /6. 0 

2 -!T3*T4+T4=»T1+TI=«‘T3)<'CJS! K#NQP)/2.0 

3 +!T4#Tl+T4*T2+Tl*T2)*CJS!K#NQ)/2,0 

4 -!Tl«T3<-T2*T3 + T2*Tl)*=CJS( K,NQM) /6.0 
OJPHI=!T2*T3 + T2<‘T4+T3=»T4)#CJPHI!K,NQ2)/6.0 

2 -!T3=>T4+T4#Tl + Tl*T3J<'CJPHl!K,NQP)/2.0 

3 ♦(T4#Tl+T4*T2+Tl*T2)*CJPHI!K,NQ)/2.0 

4 -!Tl«T3+T2*T3+T2*Tl)*CJPHI!K#NQM)/6.0 
GO TO 62 0 

CONTINUE 



TWO POINT LAGRANGIAN INTERPOLATION 

DJS = CJS!K#NQ I- CJS!K#NQM) 
DJPHI=CJPHI!K,NQ )-CJ PH I ! K , NQM I 
GO TO 620 
603 CONTINUE 
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, 











c 

c 



601 



C 

c 

c 



604 

C 

C 

c 



620 

C 

c 

c 

c 

c 

c 

600 

650 

C 

C 

C 



9610 

C 

c 

c 

500 

c 

c 

c 

c 

9630 



9624 

112 

9724 

113 



THREE POINT LAGRANGIAN INTERPOLATION 

0JS = (T4<-T3) ♦CJS(K,NQP}/2.0-(T2-«-T4)*CJS(K,NQ) 

2 +(T2<-T3I*CJS( K,NQM)/2.0 

DJPHI=(T4+T3)#CJPHI (K,NQP)/2.0-(T2<-T4J*CJPHI( K,NQ > 

2 4-CT24-T3)#CJPHI(K,NQM)/2.0 

GO TO 620 
CONTINUE 

IF (N.EQ.NQP) GG TO 600 
IF (N.EQ.NQ2) GO TO 604 

THREE POINT LAGRANGIAN INTERPOLATION 

0JS=<T2^-T3)<'CJS(K,NQ2)/2.0-(T3-i-Tl )*CJS(K,NQP) 

2 +(T24-T1)4cCJS( K,NQ)/2.0 

DJPHI=(T2+T3)*CJPHI (K,NQ2)/2. 0-( T3-«-TU*CJPHI ( K,NQPJ 
2 +IT24-T1 )«CJPHI(K,NQ)/2.0 

GO TO 620 
CONTINUE 

TWO POINT LAGRANGIAN INTERPOLATION 

DJS * CJS(K,NQPI- CJS(KfNQ) 

DJPHI=CJPHI (K,NQP)-CJPHI(K,NQ) 

CONTINUE 

CONTRIBUTION OF THE K-TH SPACE SAMPLE 

HSF(N) = HSF(N)-»-RO(K)#(SI (K)#OJS-*-OJPHI)«OS(K) 

NEXT SPACE SAMPLE 

CONTINUE 

CONTINUE 

FINAL VALUE OF THE N-TH TIME SAMPLE OF THE FIELD 

HSF(N)=C4#HSF(N) 

X(N)=N*CP 

WRITE(6,9610) N,HSF<NI 
F0RMAT(5X, I5,5X,E12.5» 

NEXT TIME SAMPLE 

CONTINUE 

END OF COMPUTATION OF THE FIELD 

PLOT BACK-SCATTERED FIELD VERSUS TIME! L IGHT-METERS) 

WRITE(6,9630) 

FORMAT(lHl) 

MOOCUR=0 

CALL PLOTP(X,HSF, ITS,MODCURJ 
IF (INC.EQ.U GO TO 112 
WRITE (6,96241 

FORMAT!// ,14X.» NORMALIZED BACK 

2 14X,»RAHP RESPONSE* ,//,14X, 

3 *X-AXIS TIME (C*OT)* ,10X,*Y-AXIS FIELD 
GO TO 113 
CONTINUE 

WRITE (6,9724) 

FORMAT!//, 14X,* NORMALIZED BACK SCATTERED FIELD*,//, 

2 14X,*IMPULSE RESPONSE*, //,14X, 

3 ‘X-AXIS TIME (C*DT)*,10X,*Y-AXIS FIELD 
CONTINUE 
RETURN 
END 



SCATTERED FIELD*,//, 
(H*C*R) * ) 



(H*R)* ) 
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3jc 4c 4: :jc ^ 3{C 3jc :(c :4c :{c :(c4c ^ ^ qjc A 

* 4c 

4c SUBROUTINE SHAPE * 

# ♦ 

4c4c 4c 4c 4c 4c 4c 4c4c 4: 4c 4c4c 4e 4c4c 4c 4c 4c 4c 4c 4c 4c 4: 

PURPOSE: 

4c 4(4: 4: 4c 4c 4s ♦ 

THIS SUBROUTINE COMPUTES THE SHAPE OF THE 
BODY FROM ITS APPROXIMATE COMPUTED PHYSICAL-OPTICS 
RESPONSE. 

IT FINDS THE CONTOUR FUNCTION RO(Z) AND DIVIDES IT 
INTO APPROXIMATELY EQUAL SEGMENTS (OS). 

USAGE: 

4c4c4c4c4c4c 

CALL SHAPE ( DS, I TER t CP i DSO , IT ) 

THE DATA NEEDED AT THE INPUT ARE DSOtAND HSF. 

PARAMETER DEFINITION: 

4c 4c 4: 4c 4c 4c 4c 4:4: 4c4c 4:4: 4c 4c 4c 4c 4c 4c 4c 4e 

ZMAX: THE MAXIMUM LENGTH OF THE BODY ON THE Z-AXIS. 
(CORRESPONDING TO THE POINT WHERE THE RAMP RESPONSE 
CHANGES ITS SIGN) 

IT: THE TRANSIT TIME NEEDED TO OBTAIN THE REFLECTION 
FROM THE POINT LOCATED AT ZMAX. (IT IS THE NUMBER 
OF TIME INCREMENTS). 

DSO: OPTIMAL LENGTH OF A SEGMENT (NEEDED FOR THE 
COMPUTED SHAPE SUBDIVISION) .(METER) 

ALL OTHER PARAMETERS ALREADY DEFINED IN MAIN PROGRAM 

REMARKS : 

IN THIS PROGRAM THE RAMP RESPONSE USED ASSUMES A 
VERTICAL H FIELD INCIDENT (HORIZONTAL POLARIZATION) 
THIS GIVES A POSITIVE INITIAL RESPONSE. 

FOR VERTICAL POLARIZATION THE INITIAL RESPONSE IS 
NEGATIVE. (THE WHOLE RESPONSE CHANGES ITS SIGN) 

SINCE THE RADIUS OF THE BODY AT EACH POINT IS 
PROPORTIONAL TO SQRT OF THE INITIAL RAMP RESPONSEt 
IT MUST BE ENSURED THAT THE RAMP USED HAS PROPER SIGN. 

METHOD OF SOLUTION: 

4c 4c 4:4c 4c 4c 4c 4c 4c 4c 4c 4c4c 4: 4c 4c 4c 4c 4c 

FIRST STEP IS TO FIND ZMAX, WHICH 
CORRESPONDS TO THE POINT WHERE THE RESPONSE 
CHANGES ITS SIGN. 

THE CONTOUR FUNCTION OF THE BODY RO(Z) IS 
PROPOkriONAL TO THE SQRT OF THE PHYSICAL-OPTICS 
RAMP RESPONSE. 

THE BODY IS DIVIDED INTO APPROXIMATELY EQUAL 
SEGMENTS, AND THE SPACE SAMPLES DELIMITING THE 
SEGMENTS ARE FOUND USING A 'CUT AND TRY METHOD' 

WHICH CONVERGES TO THE DESIRED LENGTH OF SEGMENT. 

THE RADIUS AT THE SPACE SAMPLE POINTS ARE 
INTERPOLATED FROM THE KNOWN SAMPLES OF THE RAMP 
RESPONSE. 

THE DATA OBTAINED Z , RO , OS, NS, I S USED FOR THE NEXT 
COMPUTATION OF THE FIELD. 
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c 

c 

SJBROUTINie SHAPE ( DS , ITER , CP ,□ SO » I T) 

DIMENSION OS{20) »RO(21)t Z( 2 1 ) , HSF{ SO ) 

COMMON HSF,ZtRO»ITStNS,NSl 
COMMON XMAX 
C 

C START PROGRAM 

C 

C 

C FIND LENGTH OF BODY (ZMAX ) AND ITS CORRESPONDING 

C TIME OF APPEARANCE IN THE RAMP RESPONSE (IT) 

C 

03P=DS0 

DO 710 N=l, ITS 
IF{HSF(N) .LT.0.0) GO TO 720 
710 CONTINUE 

720 CONTINUE 

IT =N-1 

AO=l./{ l.-HSF(N)/HSF{ IT )) 

ZMAX=(IT -2+AD)*CP/2. 

C 

C 

C COMPUTE SHAPE OF BODY UP TO ZMAX. 

C 

54 CONTINUE 
Z(1 )=0.0 

RO (1)=SQRT{2*HSF(2) ) 

DO 50 K=l,19 
K1 = K+1 
U=DSP/2.0 

55 CONTINUE 
NN=1 

C 

C ADJUST Z TO OBTAIN PROPER OS 

C 

Z(K1)=Z(K)+U 
53 CONTINUE 

IF(Z(K1).GT.ZMAX)G0 TO 56 
GO TO 57 

56 CONTINUE 
Z(K1)=ZMAX 

RO(K1I=0.0 
GO TO 51 

57 CONTINUE 
C 

C COMPUTE CORRESPONDING RADIUS { RO ) OF THE BODY 

C 

G=2*Z(K1) /CP+2. 

M=IFIX(G) 

M1=M+1 

RQ2=(HSF{ Ml)-HSF(M) )*( G-M)+HSF( M ) 
R0(K1)=SQRT{R02*2.) 

C 

C COMPUTE LENGTH OF SEGMENT OS 

C 

51 CONTINUE 

OS (K)=SQRT( ( Z(K1 )-Z(K) )*=cc2 + ( R0(K1)-R0(K) ) 

A=0S(K)/0SP 
C 

C ALL OBTAINED SEGMENTS MUST BE WITHIN 1/100 OF DSO 

C 

IF (A.GT.1.01) GO TO 58 
IF( Z(K1).EQ.ZMAX)G0 TO 60 
IF (A.LT.0.99) GO TO 53 
GO TO 50 
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ADJUST Z AND GO BACK TO CHECK FOR DS 



C 

c 

c 

58 



90 



50 

C 

C 

C 

60 



C 

C 



61 



30 

62 



70 

C 

C 

C 

C 

c 

9601 



9602 

40 

C 



CONTINUE 

NN=NN+1 

IF (A.GT.1.0) GO TO 
Z(KU=Z(K1)+U/NN 
GO TO 53 
CONTINUE 

Z{Kl)=Z(Kl)-U/NN 
IF (Z(Kl).GT.Z(K) ) 
Z(KI)=Z(K1)+U/NN 
NN=NN+1 
GO TO 90 
CONTINUE 

NEXT SEGMENT 



CONTINUE 

NS1=K1 

NS=NSl-i 



90 

GO TO 53 



ADJUST FOR LAST SEGMENT 

U2=0S( NS) /DSP 
IF(U2.LT.0.10) GO TO 62 
IF(U2.LT.0*95) GO TO 61 
GO TO 70 
CONTINUE 

U3=(Z(NSU-Z( NS) )/(NS-l) 

NSl=NSl-l 
NS=NSl-l 
00 80 K=l ,NS 
Kl=K+l 

Z(K1)=Z( K1)+K*U3 
Z(NS1)=ZMAX 
G=2#Z(Ki) /CP+2. 

M=IF1X(G) 

M1=M+1 

R02=(HSF(Ml )-HSF(M) )#( G-M )+HSF( M) 

R0(Kl)=SQRT(R02*2. ) 

RO(NS1)=0.0 

DS (K) = SQRT(( Z(Kl )-Z(K) ) **2 + ( RO< K 1 )-RO( K ) ) **2 ) 

CONTINUE 
GO TO 70 
CONTINUE 
Z(NS)*ZMAX 
R0(NS)=0.0 
NS2=NS-1 

0S(NS2I=SQRT( (ZMAX-Z(NS2) )«*2+R0(NS2)**2) 

NS1=NS1-1 
NS =NSl-l 
CONTINUE 

SHAPE OF BODY COMPLETED AND SUBDIVIDED 
WRITE NUMBER OF SEGMENTS OBTAINED 
WRITE (6,9601 ) ITER, NS 

FORMATdHl ,5X ,»NEW NUMBER OF SEGMENTS AFTER ITERATION* 
2 ,14,* IS* ,14, ///,10X,*Z* , IOX,*RO* ,10X,*DS* ,/////) 

DS (NS1)=99.9 
DO 40 N=l,NSl 

WRITE(6,9602) Z( N ) , R0( N) , DS (N) 

FORMAK 5X,3F10*5) 

CONTINUE 
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603 



604 



PLOT THE CONTOUR FUNCTION OF THE BODY 

NS2=NSl+l 
Z(NS2)= XMAX 
RO (NS2»= Z<MS2) 

WRITE(6, 96031 
FORMAT! IHU 
M00CUR*0 

CALL PL0TP( Zf R0tNS2fM0DCURI 
WRITE (6f9604) ITER 

FORMAT!// ,20X, 'SHAPE OF BODY AFTER ITERATION NUMBER », 
2 14, //,20X, 'X-AXIS Z Y-AXIS RADIUS !R0)») 

END OF SUBROUTINE SHAPE 

RETURN 

END 
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EXAMPLE OF DATA CARDS 



0.0 


0.0 


0.02185 


0.20791 


0.08645 


0.40674 


0.19098 


0.58779 


0.33087 


0 . 74314 


0 . 50 


0. 86603 


0.69098 


0.95106 


0.89547 


0 . 99452 


1.10453 


0.99452 


1.30902 


0.95106 


1.50 


0. 86603 


1.66913 


0.74314 


1.80902 


0.58779 


1.91355 


0.40674 


1.97815 


0.20791 


2.0 


0 . 0 
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APPENDIX F 



RELATIONSHIP BETWEEN RAMP RESPONSE AND SHAPE 



The validity of (5.2) will be demonstrated analytically 
here in the case of a metallic body of revolution. We first 
express the axially directed incident ramp field as: 



H^(z,t) 



/ t - z/C t > z/C 

^ 0 t < z/C 



(F.l) 



The physical optics current induced on the body is given by: 



J(r,t) 




X H^(r,t) 



(F.2) 



The backscattered far-zone field is given by equation (3.8) 
and its scalar equivalent for the body of revolution is 
expressed by (3.16). Substitution of the current in (.3.16) 
gives the physical optics ramp response as follows: 



H“(t) = 



4r C 
o 



/ ^[-4H^(z' ,t) sin a'lp’ds’ (F.3) 



The time derivative of the incident field expressed by (F.l) 
is unity: 



K '4 

3t at 



1 t < z/C 



(F.4) 



Substitution of (F.4) into (F.3) gives: 
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1 



(F.5) 



H^(t) 



0 



J sin a' p'ds' 



From the geometry of the body of revolution we note that; 



dp' = sin a' ds' 



Substitution in (F.5) gives the following expression: 



H|(t) 



0 



P (z) 

I P'dp' 



1 2 , , 

P (z) 



2r C 
o 



(F.6) 



Since the area function is given by: 

2 

A(z) = 7T p (z) 

the physical optics ramp response (F.6) will be related to 
A ( z ) by 



H^(t) 



2irr C 
o 



A(z) 



(F.7) 
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